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Introduction 

While  x-ray  mammography  is  the  most  reliable  method  available  at  present  for  the  detection  of  breast 
cancer  in  screening  programs,  it  still  does  not  detect  all  cancers.  A  great  deal  of  research  effort  over  the 
past  several  decades  has  been  directed  towards  the  development  of  better  and  more  effective  imaging 
systems.  These  new  systems  must  be  designed  carefully  to  ensure  they  can  produce  images  of  the  highest 
quality  possible  for  a  specified  radiation  dose  to  the  patient. 

Fourier-based  linear-systems  transfer  theory  is  often  used  to  develop  theoretical  models  of  the  signal  and 
noise  performance  of  new  system  designs.  While  it  has  been  used  successfully  in  a  number  of  new 
system  designs,  only  relatively  simple  systems  can  be  analyzed  using  this  approach.  We  are  developing 
new  theoretical  relationships  that  will  extend  the  capabilities  of  linear-systems  theory  so  that  it  can  be 
used  in  the  design  of  increasingly  complex  systems.  In  particular,  we  are  developing  new  Fourier- 
based  transfer  relationships  for  the  description  of:  a)  complex  parallel  cascades ;  and,  b)  digital 
detector  systems.  The  result  of  our  research  will  be  a  generalized  resource  “< library  ”  of  transfer 
relationships  to  be  used  by  scientists  and  engineers  to  achieve  optimal  designs  of  new  imaging  systems 
for  mammography. 


Body 

The  following  four  tasks  correspond  to  the  four  tasks  outlined  for  Year  1  in  the  Statement  of  Work  of  our 
grant  application.  A  statement  of  our  progress  on  each  task  to  date  is  summarized  in  this  section. 

Task  1.  Develop  Expertise  in  Stochastic  Point-Process  Theory 

Stochastic  point-process  theory  has  been  adopted  as  our  primary  means  of  mathematically  representing 
the  distribution  of  image  quanta  that  collectively  form  a  quantum  image.  This  includes  both  x-ray  images 
and  optical  images.  In  our  work  this  year,  we  have  also  used  this  approach  to  represent  the  spatial 
distribution  of  particular  events  that  take  place  within  an  imaging  system,  such  as  the  distribution  of 
photoelectric  events,  or  the  distribution  of  sites  where  characteristic  x  rays  are  re-absorbed.  The  use  of 
stochastic  point-process  theory  in  medical  imaging  was  recently  pioneered  by  H.H.  Barrett,  and  he 
remains  one  of  only  a  few  such  experts  in  the  world  on  this  material. 1,2,3  Drs.  Barrett  and  Myers  are 
currently  writing  a  new  book  scheduled  for  publication  in  2001  that  describes  the  fundamentals  of  this 
method.3  The  use  of  stochastic  point-process  theory  is  complex  and  requires  more  mathematical 
expertise  in  random  processes  than  most  scientists  and  engineers  who  design  imaging  systems  can 
achieve  in  a  reasonable  time.  However,  use  of  the  transfer-theory  relationships  can  generally  be 
accomplished  by  non-mathematicians,  making  the  benefits  of  the  complex  statistical  analysis  available  to 
system  designers. 

In  the  past  year,  we  have  developed  our  own  expertise  in  the  use  of  point-process  theory  to  a  level 
adequate  for  the  analysis  of  medical  imaging  systems.  Using  this  approach,  each  quantum  is  represented 
as  a  Dirac  delta  function  distributed  somewhere  within  the  image  and  a  quantum  image  is  described  as  a 
spatial  distribution  of  many  such  delta  functions.  What  we  call  “image  noise”  is  the  impact  of  statistical 
correlations  between  these  quanta.  The  mathematical  relationship  between  these  statistical  correlations 
in  the  input  x-ray  image  and  the  output  digital  image  is  described  in  terms  of  signal  and  noise  “transfer- 
theory”  relationships. 

Task  2.  Develop  Expertise  in  Moment-Generating  Function  Theory 
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An  additional  skill  needed  for  this  work  is  expertise  in  moment-generating  functions  for  this  research. 
This  has  also  be  completed  but  was  significantly  easier  than  developing  expertise  in  random  point 
processes. 

Task  3.  Develop  Concept  of  Cross  Spectral  Density  of  Parallel  Cascades 

A  generalized  approach  to  describing  transfer  of  the  noise  power  spectrum  (NPS)  through  medical 
imaging  systems  has  been  developed  over  the  past  several  years  in  which  image-forming  processes  are 
represented  in  terms  of  a  cascade  of  amplified  point  processes.  Until  recently,  this  approach  has  been 
restricted  to  serial  cascades  only.  The  development  of  a  general  expression  for  the  cross  spectral  density 
of  parallel  cascades  was  the  major  task  identified  for  the  first  year.  This  has  been  completed  using  the 
random  point-process  approach,  and  enables  the  use  of  parallel  cascades  of  image-forming  processes.  A 
manuscript  has  been  submitted  for  publication  in  Medical  Physics4  and  a  copy  is  included  as  Appendix  I 
of  this  report.  The  final  result  of  this  work  can  be  expressed  as  a  single  equation:  the  general  expression 
for  the  cross  covariance  of  amplified  point  processes  and  an  expression  for  the  cross  spectral  density  for 
wide-sense  stationary  conditions  is  given  by  Eq.  81  in  Appendix  I. 

These  results  extend  the  generalized  transfer-theory  approach  to  include  the  description  of  more  complex 
image-forming  processes  involving  parallel  cascades  of  quantum  amplification  processes.  This  parallel- 
cascade  approach  has  been  used  to  develop  a  theoretical  expression  for  noise-power  transfer  in  a  simple 
radiographic  screen  that  includes  the  effect  of  characteristic  x-ray  reabsorption.  The  result  confirms 
earlier  work  by  Metz  and  Vybomy5  who  showed  that  reabsorption  increases  image  noise  and  decreases 
the  detective  quantum  efficiency  (DQE)  at  some  spatial  frequencies.  Use  of  the  transfer-theory  approach 
facilitates  a  straightforward  generalization  to  many  new  digital  imaging  systems  including  conventional 
angiographic  and  active-matrix  flat-panel  systems. 

Task  4.  Design  and  Construct  Mono-Energetic  X-Ray  System 

Task  4  was  the  development  of  a  mono-energetic  x-ray  source  based  on  a  secondary  lanthanum  target. 
This  has  been  designed  and  construction  completed.  This  system  is  required  to  generate  experimental 
proof  of  the  correctness  of  the  results  described  in  Task  3  for  reabsorption  of  characteristic  radiation  in  a 
digital  imaging  system  using  a  Csl  input  phosphor.  During  year  n,  experiments  will  be  performed  using 
both  a  conventional  x-ray  image  intensifier  system  as  well  as  a  flat-panel  detector  using  the  General 


Figure  1.  Photograph  of  mono-energetic  x-ray 
source  with  top  removed. 


Figure  2.  Schematic  illustration  of  mono- 
energetic  x-ray  source. 
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Electric  digital  mammography  system.  Figure  1  shows  a  photograph  of  the  mono-energetic  source  and 
Fig.  2  shows  a  schematic  drawing  of  the  design.  It  will  be  used  to  measure  the  DQE  of  the  Csl-based 
systems  just  above  and  just  below  the  iodine  k  edge.  The  difference  in  the  DQE  is  due  to  the  effects  of 
reabsorption. 


The  results  described  here  are  consistent  with  the  tasks  identified  in  the  approved  grant  application.  No 
insurmountable  difficulties  are  expected  to  prevent  us  from  achieving  the  tasks  identified  for  year  n. 


Key  Research  Accomplishments 

The  key  research  accomplishments  are  itemized  here,  consistent  with  the  tasks  identified  above  plus  two 

additional  accomplishments  directly  related  to  the  current  grant: 

1 .  Develop  expertise  in  stochastic  point-process  theory. 

2.  Develop  concept  of  cross  spectral  density  of  amplified  point  processes.  This  has  been  the 
missing  link  preventing  the  use  of  parallel  cascades  in  the  design  of  complex  digital  imaging 
systems  for  mammography.  A  manuscript  describing  this  result  has  been  submitted  for 
publication  and  is  attached  as  Appendix  I. 

3.  Design  and  build  mono-energetic  x-ray  source  for  experimental  verification  of  the  parallel 
cascaded  approach  to  understand  the  effects  of  self  reabsorption  in  detectors. 

4.  Completion  and  submission  of  a  paper  describing  how  linear  systems  theory  can  be  used  to 
describe  x-ray  scatter  in  terms  of  a  scatter  operator,  attached  as  Appendix  EL6 

5.  A  review  chapter  has  been  written  in  the  past  year  and  recently  published  by  the  Society  for 
Optical  Engineering.7  This  chapter  consists  of  approximately  100  pages  and  describes  the  use  of 
linear  systems  transfer  theory  in  the  design  of  digital  x-ray  detectors.  Approximately  25%  of  this 
chapter  is  based  on  research  completed  during  the  first  year  of  the  current  grant.  This 
investigator  has  been  invited  to  give  a  course  on  this  material  during  the  SPIE  meeting  in  San 
Diego,  February  2001. 


Reportable  Outcomes 

Publications 

(students  working  under  my  direct  supervision  indicated  with  •) 

1 .  •  J.  Yao,  *T.E.  Moschandreou  and  I.A.  Cunningham.  Parallel  cascades:  New  ways  to  describe  noise 
transfer  in  medical  imaging  systems,  Medical  Physics  [submitted], 

2.  T  A  Cunningham.  *M.S.  Westmore  and  A.  Fenster.  Unified  representation  of  image  blur  and  noise 
in  linear-systems  transfer  theory  using  a  scatter  operator.  Medical  Physics  [submitted]. 

3.  LA.  Cunningham.  Applied  linear-systems  theory.  Chapter  2,  Handbook  of  Medical  Imaging,  Vol  1. 


Physics  and  Psychophysics,  Eds.  J.  Beutel,  H.L.  Kundel  and  R.L.  Van  Metter,  pp.  79-159  (The 
International  Society  for  Optical  Engineering,  Bellingham,  Washington,  2000). 

Abstracts  and  Presentations 

(students  working  under  my  direct  supervision  indicated  with  •) 

1.  Understanding  radiologic  image  quality:  From  basic  concepts  to  a  practical  tool-kit  for  scientists 
and  engineers,  American  Association  of  Physicists  in  Medicine  and  World  Congress  of  Medical 
Physics  and  Biomedical  Engineering,  Chicago,  July  2000. 

2.  •  J.  Yao,  *T.  Moschandreou  and  T.A.  Cunningham.  Cross  Covariance  of  Correlated  Point  Processes 
for  use  in  Linear-Systems  Theory,  World  Congress  of  Medical  Physics  and  Biomedical 
Engineering,  Chicago,  July  2000. 

3.  *T.  Moschandreou,  •  J.  Yao  and  T.A.  Cunningham.  Use  of  the  Cross  Covariance  in  Linear-Systems 
Theory  to  Model  the  DQE  of  Detectors  with  Fluorescence  Reabsorption,  World  Congress  of 
Medical  Physics  and  Biomedical  Engineering,  Chicago,  July  2000. 


Conclusions 

The  most  important  outcome  of  the  first  year  of  progress  has  been  development  of  the  idea  of  parallel 
cascaded  of  amplified  point  processes.  In  particular,  Eq.  (81)  of  Appendix  I  describes  the  cross  spectral 
density  which  has  been  the  missing  link  preventing  prior  use  of  this  approach.  Using  it,  linear-systems 
transfer  theory  can  be  used  to  predict  the  detective  quantum  efficiency  (DQE)  during  the  design  of 
complex  x-ray  detectors  being  developed  for  digital  mammography,  to  ensure  optimal  design  of  these 
detectors  that  will  maximize  image  quality  for  any  specified  radiation  dose  to  the  patient. 
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Appendix  I:  Parallel  Cascades 

This  appendix  consists  of  a  copy  of  a  manuscript  submitted  to  Medical  Physics.  It  describes  the  concept 
of  parallel  cascades  in  linear-systems  transfer  theory  and  is  the  primary  outcome  for  year  I. 
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Parallel  cascades:  New  ways  to  describe  noise  transfer  in  medical  imaging  systems 

J.  Yao,  T.E.  Moschandreou 

Imaging  Research  Laboratories,  The  John  P.  Robarts  Research  Institute , 

London,  Ontario  N6A  5K8,  Canada 


and  I.A.  Cunningham 

Imaging  Research  Laboratories,  The  John  P.  Robarts  Research  Institute  and  Department  of  Diagnostic  Radiology, 
London  Health  Sciences  Center  and  The  University  of  Western  Ontario, 

London,  Ontario,  Canada 
(November  16,  2000) 

A  generalized  approach  to  describing  transfer  of  the  noise  power  spectrum  (NPS)  through 
medical  imaging  systems  has  been  developed  over  the  past  several  years  in  which  image-forming 
processes  axe  represented  in  terms  of  a  cascade  of  amplified  point  processes.  Until  recently,  this 
approach  has  been  restricted  to  serial  cascades  only.  Here,  we  develop  a  generalized  expression 
for  the  cross  covariance  of  amplified  point  processes  and  an  expression  for  the  cross  spectral 
density  for  wide-sense  stationary  conditions.  These  results  extend  the  generalized  transfer- 
theory  approach  to  include  the  description  of  more  complex  image-forming  processes  involving 
parallel  cascades  of  quantum  amplification  processes. 

This  parallel-cascade  approach  is  used  to  develop  a  theoretical  expression  for  noise-power 
transfer  in  a  simple  radiographic  screen  that  includes  the  effect  of  characteristic  x-ray  reabsorp¬ 
tion.  The  result  confirms  earlier  work  by  Metz  and  Vyborny  who  showed  that  reabsorption 
increases  image  noise  and  decreases  the  detective  quantum  efficiency  (DQE)  at  low  spatial 
frequencies.  Use  of  the  transfer-theory  approach  facilitates  a  straightforward  generalization 
to  many  new  digital  imaging  systems  including  conventional  angiographic  and  active-matrix 
flat-panel  systems. 

Key  words:  random  point  processes,  amplified  point  processes,  cross  covariance,  cross  spec¬ 
tral  density,  noise  power  spectrum,  transfer  theory 


Submitted  to  Medical  Physics  May  2000 

I.  INTRODUCTION 

Medical  x-ray  imaging  systems  must  be  designed  to 
ensure  that  maximum  image  quality  is  obtained  for  a 
specified  radiation  dose  to  the  patient.  While  there  are 
many  aspects  of  “image  quality” ,  one  important  con¬ 
sideration  is  image  noise  as  described  by  the  Wiener 
spectrum,  or  noise  power  spectrum  (NPS)  [1-3].  The 
NPS  describes  the  spectral  decomposition  of  second- 
moment  statistics  in  terms  of  spatial  frequencies  un¬ 
der  wide-sense  stationary  (WSS)  conditions  [2,3].  It  is 
required  for  the  determination  of  other  image-quality 
metrics  used  to  quantify  image  quality  and  system 
performance  including  the  noise-equivalent  number  of 
quanta  (NEQ)  [1,4,5],  which  describes  an  equivalent 
number  of  quanta  forming  an  image,  and  the  detective 
quantum  efficiency  (DQE)  [4-9],  which  describes  the 
ability  of  an  imaging  system  to  make  efficient  use  of 
the  incident  image  quanta. 

Over  the  past  several  years,  a  generalized  transfer- 
theory  approach  [10-12]  has  been  developed  to  de¬ 
scribe  how  the  NPS  is  transferred  from  the  input  of 
an  imaging  system  to  the  output  image  [13-18].  Of 
particular  significance  in  this  development  was  a  de¬ 
scription  of  how  the  NPS  is  transferred  through  quan¬ 
tum  gain  and  quantum  scattering  stages  by  Rabbani, 
Shaw  and  Van  Metter  [13]  and  by  Barrett,  Wagner  and 


Myers  [19,20].  This  generalized  description  of  image 
noise  has  resulted  in  a  comprehensive  frame-work  for 
the  understanding  of  system  performance  built  upon 
a  communication  theory  based  approach. 

Using  this  approach,  many  imaging  systems  can  be 
represented  in  terms  of  serial  cascades  of  three  ele¬ 
mentary  processes:  i)  quantum  amplification;  ii)  quan¬ 
tum  scattering;  and,  iii)  linear  filters.  Transfer  of  sig¬ 
nal  and  noise  through  these  models  can  be  described 
by  cascading  transfer  relationships  of  each  elementary 
process.  In  Appendix  A,  transfer  properties  of  these 
three  elementary  processes  are  summarized.  This  ap¬ 
proach  has  been  used  recently  to  describe  signal  and 
noise  transfer  and  the  DQE  of  a  number  of  x-ray 
medical  imaging  systems,  including  film  screen  sys¬ 
tems  [21-23],  active-matrix  flat-panel  systems  for  digi¬ 
tal  radiography  [24-26],  video-based  systems  for  portal 
imaging  and  radiation  therapy  verification  [27,28],  and 
other  new  system  designs  [29-32]. 

We  are  developing  a  number  of  new  transfer-theory 
relationships  to  describe  noise  transfer  through  pro¬ 
cesses  where  the  three  elementary  processes  described 
above  are  inadequate.  These  new  relationships  form 
the  basis  of  new  theoretical  “tools”  that  can  be  used 
by  scientists  and  engineers  developing  or  assessing  new 
system  designs.  Of  particular  practical  importance  is 
the  spatial- frequency-dependent  form  of  these  tools  for 
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WSS  conditions.  They  can  be  used  to  make  a  theo¬ 
retical  prediction  of  the  NEQ  or  DQE  of  a  particular 
system  design. 

One  such  current  limitation  of  the  transfer-theory 
approach  is  that  it  has  been  restricted  to  serial  cas¬ 
cades  of  the  elementary  processes.  This  excludes  the 
situation  where  more  than  one  image-forming  process 
must  be  summed  to  form  an  image  [16].  For  example, 
x  rays  interacting  in  a  radiographic  screen  generally 
do  so  by  the  photo-electric  interaction.  This  process 
often  results  in  the  emission  of  a  characteristic  x  ray 
that  may  be  reabsorbed  elsewhere  in  the  screen.  Light 
is  generated  at  both  the  primary-interaction  and  re¬ 
absorption  sites,  but  with  different  intensities.  In  ad¬ 
dition,  the  reabsorption  site  is  randomly  located  but 
spatially  correlated  with  the  primary-interaction  site. 
Light  from  both  sites  contribute  to  production  of  the 
final  radiographic  image  recorded  on  film.  However,  it 
is  not  possible  to  describe  image  noise  as  the  sum  of 
these  correlated  image-forming  processes  using  a  sim¬ 
ple  serial  cascade,  and  hence  the  effect  of  reabsorption 
has  not  been  included  in  any  transfer-theory  analysis. 

In  this  article,  we  extend  the  capabilities  of  the 
transfer-theory  approach  so  that  more  complex  sys¬ 
tems  requiring  both  serial  and  parallel  cascades  of 
these  elementary  processes  can  be  represented.  This 
is  accomplished  by  developing  a  general  expression  for 
the  cross  covariance  and  cross  spectral  density  of  noise 
processes  that  can  be  incorporated  into  the  transfer- 
theory  analysis  (see  Appendix  B).  Use  is  made  of  ran¬ 
dom  point  process  theory,  where  a  quantum  image  is 
represented  as  a  two-dimensional  spatial  point  process 
in  which  each  quantum  is  represented  as  a  point  im¬ 
pulse  [19,20,33].  A  general  expression  is  derived  for 
the  cross  covariance  of  two  correlated  random  point 
processes.  This  is  then  simplified  for  the  special  case 
of  wide-sense  stationary  random  point  processes  where 
the  cross  spectral  density  function  is  derived. 

We  derive  here  a  general  expression  for  the  cross 
covariance  of  a  pair  of  correlated  point  processes.  Of 
practical  importance  for  applications  in  medical  imag¬ 
ing  is  the  special  case  where  each  point  process  repre¬ 
sents  a  subset  of  a  common  input  point  distribution. 
The  cross  covariance  of  these  two  subsets  is  derived, 
and  then  generalized  to  describe  the  cross  covariance 
of  the  two  distributions  after  they  subsequently  un¬ 
dergo  an  arbitrary  cascade  of  quantum  amplification 
and  scattering  processes.  It  is  shown  that  a  very  sim¬ 
ple  closed-form  expression  for  the  cross  covariance  and 
cross  spectral  density  exists  under  WSS  conditions, 
where  the  cross  spectral  density  is  the  Fourier  trans¬ 
form  of  the  cross  covariance.  Use  of  the  cross  spectral 
density  is  then  demonstrated  in  an  analysis  of  noise  in 
a  radiographic  screen  with  reabsorption.  This  prob¬ 
lem  was  first  solved  by  Metz  and  Vyborny  [34]  using 
a  very  different  type  of  statistical  analysis.  Our  work 
confirms  their  result,  and  is  of  a  more  general  nature 
that  is  readily  extended  to  describe  reabsorption  in 


other  imaging  systems  including  digital  flat-panel  ra¬ 
diographic  systems. 

Throughout  the  following  description,  we  use  a  no¬ 
tation  where  the  overhead  tilde  (eg.  N)  indicates  a 
random  variable,  overline  (e.g.  q)  indicates  a  mean 
value  and  bold  face  (eg.  r)  indicates  a  vector. 

II.  THEORY 

The  framework  of  this  analysis  is  based  in  part  on 
earlier  works  by  Barrett  et  al  [19,20]  who  developed 
the  use  of  random  point  process  theory  for  studying 
noise  in  imaging  systems.  A  random  point  process  is 
any  random  process  for  which  all  sample  functions  can 
be  represented  as  a  distribution  of  points,  and  we  will 
represent  each  point  as  a  spatial  Dirac  S  function.  For 
instance,  a  quantum  image  is  described  as  a  spatial 
distribution  of  S  functions.  However,  these  points  may 
also  represent  a  spatial  distribution  of  certain  events, 
such  as  a  distribution  of  photo-electric  events,  or  a 
distribution  of  photo-electric  events  when  a  K  x  ray  is 
reabsorbed. 

A  random  point  process  is  associated  with  two  im¬ 
portant  quantities:  the  location  of  points  in  space 
where  events  occur  and  the  number  of  such  points  [33]. 
The  mathematical  realization  of  a  spatial  point  process 
can  be  expressed  as  a  sequence  of  random  impulses, 
given  by 

N 

q(r)  "  ~  rn),  (!) 

n= 1 

where  r  is  a  multidimensional  spatial  coordinate  vec¬ 
tor  in  space  S  where  the  point  process  is  defined,  rn  is 
a  continuous  random  vector  describing  the  location  of 
the  nth  point  falling  in  S  and  N  is  a  random  variable 
describing  the  number  of  points.  The  ensemble  of  ran¬ 
dom  vectors  describing  the  positions  of  all  N  points 
is  {r n  :  n  =  1, 2,  •  ■  • ,  N}.  In  this  section,  we  derive  a 
general  expression  for  the  cross  covariance  of  two  cor¬ 
related  point  processes  drawing  on  previous  work  by 
Barrett  et  al.  [19,20]. 

A  quantum  image  is  represented  as  a  sample  q(r)  of 
the  random  point  process  given  by  Eq.  (1),  where  the 
space  S  of  points  denotes  the  two-dimensional  image 
area.  Our  analysis  is  also  applicable  to  higher  dimen¬ 
sional  space.  Although  the  size  of  S  is  arbitrary,  an 
infinite  size  is  required  for  the  analysis  under  WSS 
conditions,  and  hence  we  consider  S  to  be  infinite  in 
size. 


A.  Cross  covariance  of  point  processes 

For  the  general  case,  we  consider  two  random  spatial 
point  processes, 
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&(»*)  =  j>(r-riJ),  (2) 

n—  1 

and 

nb 

&(r)  =  22f(r-rf),  (3) 

j= i 

which  may  or  may  not  be  statistically  correlated.  The 
cross  correlation  of  gU(r)  and  qB(r)  is  the  mean  of  the 
product  qA(r)q%(r')  [2],  i.e., 

R-AB(r,r')  =  (qA(r)q*B(r')),  (4) 

where  *  denotes  a  complex  conjugate  and  ( )  repre¬ 
sents  an  expectation  operator.  The  cross  covariance  of 
qA(r)  and  qB{r)  is  given  by 

K AB{r,r')  =  RAB(r,r')  -  (qA(r)} (q*B(r')) ■  (5) 

Barrett  et  al.  [19,20]  have  shown  that  the  mean  of 
qA(r)  in  Eq.  (2)  is  given  by 

qA(r)  =  [qA(r))  = /^2prf.A(r\NA)\  ,  (6) 

7  \n- 1  ”  /  Na 

where  pr^A(r\NA)  is  the  conditional  probability  den- 

sity  function  of  the  nth  point  of  the  process  qA(r)  eval¬ 
uated  at  rn  =  r  for  a  specified  value  of  NA,  and  we 
denote  by  { the  average  over  NA.  Similarly, 

<ZbM  =  (<?£f(i"))  =  (y2prfB  (r\NB)\  .  (7) 

\i= i  •  nb 

If  qA(r)  is  statistically  independent  of  qB(r ),  the  cross 
correlation  R AB(r,r')  is  equal  to  the  product  of  their 
means,  and  the  cross  covariance  K ab^,'*'*)  in  Eq.  (5) 
becomes  zero.  Spatial  point  processes  qA{r)  and  Qb{t) 
are  then  called  uncorrelated. 

In  order  to  calculate  the  mean  of  the  product 
QA(r)q*B(r')  in  Eq.  (4)  where  qA(r)  and  qB(r)  may  be 
statistically  correlated,  we  must  average  over  all  ran¬ 
dom  quantities  {fA},  {rf },  NA  and  NB  in  processes 
qA(r)  and  £s(r').  The  procedure  is  divided  into  two 
steps  [19].  The  first  one  is  to  take  the  conditional  ex¬ 
pectation  of  the  continuous  random  quantities  {fA} 
and  {rf }  for  fixed  NA  and  NB.  By  definition  of  the 
statistical  average  over  a  continuous  random  variable, 
shown  in  Appendix  C,  we  have 

E{g^(r)^(r')|iVA,iVB} 

=  f  drA- •  ■  f 

Joo  Joo  tn=1  j-l 


where  pr({r£},  {rf  }\NA,  NB)  is  the  conditional  joint 
density  function  of  random  variables  {rjJ}  and  {rf} 
given  Na  and  NB ,  and  we  use  the  symbol  to  de¬ 
note  a  multidimensional  integral  over  all  <S.  Using  the 
property  of  marginal  densities,  shown  in  Appendix  C, 
we  obtain 

Air 

=  EE  f^fdrf{s(r-ri)S(r'-rf) 

n=lj=1Jo°  Joo  l 

xprfAfB(rA,rf\NA,NB^.  (9) 

It  follows  from  the  sifting  property  of  delta  functions 
that 

=  Y.!L»if‘„t;{TyiNA’NB)'  (10) 

n=l  j= 1 

where  pr^x  f,B(r,r'\NA,NB)  is  the  conditional  joint 

*  n  »*  j 

density  function  of  fA  and  rf  for  fixed  NA  and  NB , 
evaluated  at  rA  —  r  and  rf  =  r*.  Next,  by  averag¬ 
ing  Eq.  (10)  over  NA  and  NB ,  we  obtain  the  cross 
correlation  of  qA(r )  and  qB(r)  given  by 

RAB(r,r#)  =  (ee  pr^n  ,rf{r’r'\NA’NB)) 

\n=l  j=l  'NA,NB 

(ii) 

Prom  Eqs.  (5)-(7)  and  (11),  therefore,  the  cross  co- 
variance  of  qA (r)  and  qB  (r)  is  given  by 

Kab(t,  r')  =  /^  E  P  rf?  ,rf  (r>  r'\NA>  nB )\ 

\n=lj=l  f^A^B 

-  (j2vTrt{r\NA)\  l'^2pirf{r'\NB)f  ■  (12) 

\n=  1  ”  / Na  \  j=1  '  Na 

Without  loss  of  generality  we  assume  that  each  point 
has  the  same  conditional  probability  density  function 
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for  given  NA  or  NB .  The  probability  density  functions 
are  therefore  independent  of  the  indices  n  and  j  and 
we  simplify  our  notation  by  using  rA  and  rB  instead 
of  rA  and  rf.  Thus,  Eq.  (12)  becomes 

KAB(r,r')  =  (NANBpii.A  fB(r,r'|lVA,  ATB)\  _ 

\  *  V  '/NA,NB 

-^prfA(r|^)^JVBprfB(r'|Ars)^  (13) 

where  the  NA  x  NB  terms  in  the  double  sum  over  n 
and  j  are  identical. 

Equation  (13)  is  a  general  expression  for  the  cross 
covariance  of  two  random  point  processes.  It  is  valid 
for  both  stationary  and  non-stationary  random  pro¬ 
cesses.  For  imaging  applications  where,  in  general, 
NA,NB  >  1  and  the  probability  density  and  the  joint 
density  functions  are  independent  of  NA  and  NB ,  the 
statistical  nature  of  NA  and  NB  can  often  be  ignored 
and  the  cross  covariance  of  two  quantum  images  is  then 
given  by 

KAB(r,r')&NANB  |pr f  B(r,  r')  -  prf  x(r)prfB(r')] . 

(14) 

B.  Cross  covariance  of  random  subsets  of  a 
random  point  process 

If  the  quanta  in  two  images  are  independent  of  each 
other,  the  cross  covariance  of  the  two  images  will  be 
zero.  This  is  certainly  the  case  when  two  images  are 
acquired  independently  of  each  other.  However,  we  are 
interested  in  the  special  case  where  two  point  distri¬ 
butions  (images)  are  not  independent.  A  simple  ex¬ 
ample  of  this  occurs  when  the  two  point  distributions 
are  each  random  subsets  of  a  common  input  point  dis¬ 
tribution  or  image.  If  quanta  in  the  input  image  are 
statistically  correlated,  there  will  in  general  be  a  non¬ 
zero  cross  covariance  between  the  two  subsets.  In  this 
section,  the  cross  covariance  of  two  point  distributions 
is  determined  when  each  represents  a  random  subset 
of  a  correlated  input  image. 

The  process  of  randomly  selecting  points  from  a  dis¬ 
tribution  is  illustrated  in  Fig.  1.  This  random  point- 
selection  process  represents  a  sequence  of  independent 
trials  in  which  each  trial  makes  a  random  determina¬ 
tion  for  each  point  in  the  input  distribution.  The  point 
is  selected  to  path  A  with  probability  £,  and  path  B 
with  probability  £.  That  is,  each  trial  is  described  in 
terms  of  two  binomial  random  variables,  denoted  by 
£n  and  (n  for  the  nth  trial,  where  each  random  vari¬ 
able  can  have  a  value  of  0  or  1  only.  Each  trial  is 
independent  of  all  others,  but  we  will  allow  statistical 
relationships  between  variables  |n  and  for  a  given 


u(r) 


FIG.  1.  Illustration  of  randomly  selecting  points  from  a 
random  point  process. 

trial  n.  For  N  trials,  the  set  of  random  variables  £n 
and  £n  required  to  describe  a  single  image  is  given  by 
{in,Cn--n  =  l,2,---,N}. 

Similar  to  Eq.  (1),  the  input  is  a  general  spatial  point 
process  rewritten  as 

N 

qin(r)  =  '£/6(r-Rn),  (15) 

n—1 

where  Rn  is  a  random  vector  describing  the  position 
of  the  nth  quantum  in  the  input  image.  The  point  se¬ 
lection  process  is  assumed  to  be  independent  of  qin(r)- 
The  outputs  <£4  (r)  and  <2s(r)  can  therefore  be  written 
in  terms  of  the  random  variables  and  Cn  as 

'  N  _ 

U (r)  =  ^2£nd(r  -  Rn) 

<  "l1  •  (16) 

N 

Qb{t)  =  53  C nS{r  -  Rn) 

<  n=  1 

In  Eq.  (16),  the  point  processes  gUM  and  Qb( r)  are 
random  subsets  of  ^n(r).  We  are  interested  in  the 
cross  covariance  of  qA{r)  and  ^(r).  In  the  following, 
we  examine  both  first-order  and  second-order  statistics 
of  qA{r)  and  qsir). 


1.  Mean 


We  calculate  the  mean  of  gUM  m  Eq.  (16),  in  two 
steps.  The  first  step  is  to  take  the  conditional  expec¬ 
tation  of  qA{r)  for  fixed  {Hn}  and  AT,  i.e.,  for  fixed 
qin{ r)  refering  to  the  sample  space  of  qin(r)i  and  then 
to  average  over  qin  (r).  Thus,  we  have 
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/  N 

E{<jU(r)km(r)}  =  /  ^  £nS(r  -  Rn) 


n=l 
N 


®»(r) 


^E{|„}<5(r-Rn) 


n=l 

N 


Qin  (T) 


=‘(^(r-Rn)  =  iqin(r).  (17) 


n=l 


where  |n  has  the  same  mean  value  for  all  n,  given  by  C> 
and  we  use  the  symbol  |9in(r)  to  denote  the  computa¬ 
tion  condition,  i.e.,  for  fixed  &n(r).  Next,  by  averaging 
over  qin(r)i  we  can  obtain  the  mean  of  qU(r)  given  by 


points  in  A  and  B,  then  Cn  is  orthogonal  with  fn,  i.e., 
Cn  =  (1  -  |n)  and  =  0.  With  this  notation,  we 
have 


E{^(r)^(r/)|gin(r)}| 

N 

=  RCC  ^  £(r  -  Rn)S(r'  -  Rn) 


n=l 


m=i 


—  ■R'CC  qiniy )qin(y  ) 

I  n—j 

When  n  ^  j,  random  variable  |n  is  independent 
and  E{CnCj}  =  C  C*  In  this  case,  we  have 


(22) 


of  0 


9j4(r)  =  E{^(r)}  =  f  ®»(r). 

Similarly, 

gB(r)  =  E|gs(r)  j  =  CQin(r). 


(18) 


E{^(r)«B(0|?in(r)} 

/  AT  N  _  _  _  „  \ 

(19)  =  (  E  E  -  ~  ^i)  Qin(r)\  _ 

'n=l  7  =  1 


2.  Cross  correlation  and  cross  covariance 


n=l  j=l 


The  cross  correlation  of  qA(?)  and  2b  (Oj 
RAB(r,r'),  is  the  mean  of  the  product  (**')> 

i.e., 

R>iB(r,T-')  =  e|  qA(r)q*B(r')} 

=  (E  -  R„)  E  C i  &(r'  -  Rj))-  (20) 

\n= 1  j=l  / 


=  tCqin(r)qin(r') 


n^j 


(23) 


Adding  Eqs.  (22)  and  (23),  the  conditional  expectation 
of  the  product  qA(r)q*B(r')  is  given  by 


E|gA(r)9B(r')ktn(r)} 

=  R?c  qin(r)qin(r')\n_.  + £  CQin(r)qin(r')\n^  (24) 


Again,  using  an  approach  similar  to  Barrett  et  al  [19], 
computation  of  the  expectation  in  Eq.  (20)  is  divided 
into  two  steps.  That  is,  the  first  one  is  to  average  over 
{|n}  and  {Cn}  for  fixed  input  qin{r ),  and  then  average 
over  QinW*  There  are  two  cases  to  be  considered  in 
Eq.  (20),  corresponding  to  n  =  j  and  n  ±  j.  When 
n  =  j,  which  has  N  terms, 


E{gA(r)^(r,)km(r)| 

^  J  n—j 

"  „  _  \ 

E)  inC n  3(r  ~  Rn)$(r  ~  Rn)  qin(r)  / 

n=l  '<«.< 


•  (21) 


n=j 


which  after  averaging  over  (ji„(r)  yields 
Rab(**,t-')  =  R^Rm(r,r')|  .  +  £CRin(r,r') 

ln=j 


n^j 

(25) 


where  Ri„(r,r-')  =  E{?in(r)^(r')}  is  the  autocorre¬ 
lation  of  qin(r)  and  is  given  by 


Rin(r.r')  =  Rin  0",r')  +Rin(r,r') 


n=3 


n^j 


(26) 


Equation  (25)  can  be  further  simplified.  Note  that  for 
a  general  point  process  as  in  Eq.  (15),  we  have  [19,20] 


It  is  convenient  to  denote  the  cross  correlation  of  Cn 
and  Cn  f°r  n  —  1,2, •■•,1V  as  =  E{CnCn}-  It 
is  non-zero  when  the  two  random  variables  are  non- 
orthogonal.  For  instance,  if  the  two  images  A  and 
B  represent  identical  subsets  of  the  input  distribution, 
where  £n  =  |n>  they  are  correlated  and  the  cross  corre¬ 
lation  of  Cn  and  Cn  is  R^.  If  the  two  images  represent 
complementary  subsets,  where  there  are  no  common 


Rin(r,r')|  =qin(r)5(r-r').  (27) 

\n=j 

From  Eqs.  (25)— (27),  therefore,  the  cross  correlation 
of  qA  (»*)  and  qs  (r)  becomes 

R^B(r,r')  =  (Rcc  -  {ft qin(r)S(r  -  r')  +  |CR in{r,r') 
=  KiC  qin(r)S(r  -  r')  +  £<Ri„(r,r'),  (28) 
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where  is  the  cross  covariance  of  random  variables 

in  and  Cn,  given  by  =  R*<  -  £  C 
Finally,  the  cross  covariance  of  ^(r)  and  ^(r)  is 
given  by 

K AB{r,r')  =  R,u?(r,r')  -  qA(r)qB{rf) 

=  K €C«n(rM(r  -rO  +  fCK4n(r>r#),  (29) 

where  Kin(r,r')  is  the  autocovariance  of  the  input 
point  process  <?i„(r). 

Equation  (29)  is  a  general  expression  for  the  cross 
covariance  of  ^(r)  and  <?£(r),  where  each  is  a  random 
subset  of  the  input  point  process.  It  consists  of  two 
components.  The  first  represents  uncorrelated  noise 
power  given  as  a  <5  function  scaled  by  the  cross  covari¬ 
ance  of  the  two  binomial  random  variables  and  the 
mean  number  qin(r)  of  quanta  per  unit  area  in  the  in¬ 
put.  This  component  is  zero  when  A  and  B  represent 
independent  subsets  of  the  input,  and  non-zero  oth¬ 
erwise.  The  second  component  represents  correlated 
noise  and  is  proportional  to  the  cross  covariance  of  the 
input  point  process  Kjn(r,r'). 


C.  Cross  covariance  following  an  amplified  point 
process 

A  more  general  case  involves  the  cross  covariance  of 
two  point  processes  that  undergo  an  amplified  point 
process  subsequent  to  selection  as  illustrated  in  Fig.  2. 
Following  the  work  of  Rabbani,  Shaw  and  Van  Metter 
[13],  and  Barrett  et  al  [19,20],  an  amplified  point  pro¬ 
cess  is  considered  to  be  a  random  point  process  where 
each  point  is  converted  into  a  random  “cluster”  of  kn 
secondary  points  randomly  distributed  by  the  random 
vectors  { Ank  :  k  =  1, 2,  •  •  • ,  kn},  which  is  the  cascade 
of  the  elementary  processes,  quantum  gain  and  quan¬ 
tum  scatter,  described  in  Appendix  A.  If  each  ampli¬ 
fication  process  is  independent  of  all  others,  the  point 
processes  qU(r)  and  qn  (r)  in  Fig.  2  can  be  written  as 

QA(r)  =  E&£*(r  -  -  ^nfc) 

<  nZX  *=;  ,  (30) 

?BM  =  £C«E^-Kn-A*) 

w  n=l  k= 1 

where  the  nth  input  quantum  randomly  located  at  Rn , 
if  passed  to  paths  A  and  B  produce  k£  and  k%  quanta 
respectively,  and 


nk 


=Rn  + 


and 

A  o  i  a  ^ 
^ nk — d-  ^ nk 


t APP:  Amplified  Point  Process 

FIG.  2.  Two  point  processes  following  an  amplified 
point  process  subsequent  to  selection  of  a  random  point 
process. 

are  random  vectors  describing  the  positions  of  the  kth 
quantum  produced  by  the  nth  primary  in  the  image 
plane  for  each  path.  As  assumed  by  Barrett  et  al. 

[19,20],  we  assume  the  random  displacements  Anfc  and 
~  B 

Ank  are  independent  of  all  others,  possibly  depending 
on  the  position  of  the  primary  for  non-stationary  pro¬ 
cesses. 

Denote  by  prd(A^fc|i?„)  the  univariate  probability 
density  function  of 

^  An£  :  k  =  1, 2,  •  •  • ,  kn]  n  =  1, 2,  •  •  • ,  n) 

given  the  primary  Rn.  The  following  relationship  is 
then  known  [19,20], 

pr^AjJjJln)  =  [*»(*»)]  (Anfc  +  Rn,  Rn)  , 

(31) 

where  the  mean  number  k£{Rn)  of  secondaries  result¬ 
ing  from  primary  Rn  in  path  A  is 

kn(Rn)  =  f  drp$(r,Rn),  (32) 

J  OO 

and  p£(r,  Rn),  for  path  A,  is  defined  as  the  mean  dis¬ 
tribution  of  secondaries  at  r  when  a  primary  is  ab¬ 
sorbed  at  JR„.  Similarly  for  path  B,  we  have 

prf(A®  |Hn)  =  [^(i?„)]~V£(A**+i*n,i*n), 

(33) 

where  pr  I(A£*|J*n  )  is  the  univariate  probability  den¬ 
sity  function  of  {Anfc}  given  jRn,  and 

kn(Rn)  —  (  drp%(r,R„).  (34) 

J  OO 

We  now  examine  the  means  and  cross  covariance  of 
qA(r )  and  <7s(r)  for  these  amplification  processes. 
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1.  Mean 

The  mean  of  the  point  processes  g^(r)  and  qB(r)  are 
calculated  using  an  approach  similar  to  that  described 
by  Barrett  et  al.  [19].  For  path  A,  the  procedure  is  as 
follows. 

~  A 

(a)  Average  over  displacements  {AnAJ  for  fixed  Rn ; 


E(C){^W}  =  E(c){E(6){9A(r)}|  =  f2tpA(r,Rn). 

{  1  n—l 

(38) 

Step  (d),  for  the  continuous  random  variables  {Rn} 
with  the  univariate  conditional  density  function 
pr^  (jRn|iV)  given  N ,  requires  the  computation  of 


(b)  average  over  number  k£  of  secondaries  for  fixed 

Rn j 

(c)  average  over  binomial  random  variables  £n  for 
fixed  qin(r ); 

(d)  average  over  positions  {Rn}  of  input  quanta  for 
fixed  N ;  and, 

(e)  average  over  total  number  N  of  input  quanta  to 
get  the  result  E{^(r)}. 

Steps  (a)  and  (b)  calculate  the  conditional  expectation 
of  secondary  points  given  a  primary.  Step  (c)  averages 
over  the  selection  of  quanta  entering  path  A  for  fixed 
input.  An  average  over  the  input  point  process  is  ob¬ 
tained  in  steps  (d)  and  (e).  In  an  attempt  to  simplify 
the  notation,  we  will  use  the  step  label  as  a  subscript 
to  express  the  result  of  a  step  (a)-(e).  For  example, 
denote  by  E(a){gU(r)}  the  result  of  step  (a). 

Step  (a)  is  the  statistical  average  over  the  continu¬ 
ous  random  vectors  { A^fc}  given  Rn.  In  a  similar  way 
to  Eqs.  (8)— (10),  we  can  obtain 

N 

E(tt){flU(r)}  =  £*»£[*„(*»)]  Pa(r,Rn),  (35) 

n—l  k-1 1 


E(d){?A(r)} 

—  f  dR\  *  *  •  f  dRjsj  ^  £  Pd  (t >  Rn) 

Joo  Joe  n=1 

=  2  £  [  dRn  Pdfr ,  Rn)  P*Rn{Rn\N) 

n= 1  J(X> 

N  r 

=  £</ 
n=l  Joc 


dRPd(r,R)prklR\N), 


(39) 


where  the  property  of  marginal  densities  was  invoked 
and  the  integration  variable  Rn  is  renamed  R.  Since 
the  conditional  expectation  of  qin(r)  for  fixed  N  is 
given  by  (see  Eq.  (6)) 


E{9m(r)|Ar}  =  ^  prkn(r\N),  (40) 

n= 1 


then  Eq.  (39)  becomes 

E(d){^(r)}  dRpj(r,R)  E{qin(R)\N}. 

(41) 


where  the  probability  density  function  pr^(A^fc|iln) 
has  been  expressed  in  terms  of  p^{r^Rn)  from 
Eq.  (31).  We  can  simplify  Eq.  (35)  into 

N  _l 

Bw{Mr)}  =  i;^[«(iln)]  Pd(r,Rn),  (36) 

n—\ 

since  the  kernel  in  the  sum  over  k  is  independent  of  k. 

Step  (b)  requires  the  average  of  Eq.  (36)  over  dis¬ 
crete  value  k£  given  Rn.  This  leads  to  cancellation 
of  [k£ (Hn)]_1  in  Eq.  (36)  because  of  the  conditional 
expectation  of  k£  given  Rn  as  shown  by  E(6){fc£}  = 
k£(Rn).  Thus 

E(6){?A(r)}  =E(6)|E(0){^(r)}}  =  Y^{r,Rn). 

(37) 

Step  (c)  can  be  obtained  simply  by  replacing  £„  in 
Eq.  (37)  with  i.e., 


Step  (e),  the  average  of  Eq.  (41)  over  N,  yields 

E{fo(r)}  =ij  dRpd(r,R)  qin(R),  (42) 

where  E{£jn(.R)|iV}  is  averaged  over  N.  Since 

qAo(r)  =  |ftn(rj,  where  point  process  <jUo(r)  is  the 
output  of  the  point-selection  process  for  path  A  (see 
Fig.  2),  we  obtain 

K{qA(r)}=  f  dRpj{r,R)qAo{R).  (43) 

*  *  Joo 

Equations  (42)  and  (43)  are  general  expressions  for 
the  mean  of  an  amplified  point  process.  Similarly,  the 
mean  of  <?#(r)  is  given  by 


v{<iB(r)}=cf  dRp%(r,R)qin(R ) 

=  f  dRpd  (r ,  R)  qBo  (R)-  (44) 

J  OO 
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2.  Cross  Correlation  and  Cross  Covariance 


We  now  calculate  the  cross  correlation  for  the  out¬ 
put  point  processes  qA(r)  and  Qb{^)-  By  definition, 
R,4B(r,r')  is  given  by 


RAB(r,r')  =  E{qA(r)q*B(r')} 

j  N  A  N  \ 

\n=l  k= 1  n=l  k=  1  / 


(45) 


N 


n=l 


(c) 


N 

=  YjE{^n)  Pd(r>Rn)Pd(r'iRn) 
n—  1 

N 

=  R€<E34(r,-Rn)p*(r,,'Rn)- 

n=l 


(50) 


The  calculation  for  Steps  (d)  and  (e)  on  Eq.  (50)  is  now 
similar  to  what  was  done  in  Eqs.  (39)-(42).  Thus,  the 
result  is  given  by 


Similar  to  the  computation  of  the  mean  of  qU(r),  we 
calculate  the  expectation  in  Eq.  (45)  by  the  five  steps 

(a)-(e)  shown  above.  Step  (a)  is  to  average  over  dis- 
~  A  ~  B 

placements  {Anfc}  in  qU(r)  and  {Anfc}  in  gfi(r)  for 
fixed  {Rn},  denoted  by  E(a){^(r)qfB(^')}-  Step  (b) 
is  to  average  over  k£  and  fc®  for  fixed  {fin},  denoted 
by  E (b){QA(r)qB(r')}-  We  assume  the  point  ampli¬ 
fication  processes  in  paths  A  and  B  may  depend  on 
incident  locations  {fin},  but  are  independent  of  all 
other  terms.  That  is,  both  the  gain  factors  {k^}  and 
{k%}  are  independent  for  all  n  and  the  scatter  vectors 

{Anfc}  and  {Anfc}  are  independent  for  all  n  and  fc. 
Therefore,  we  can  write  down  the  results  for  steps  (a) 
and  (b)  as 

E(a){4U(r)<j'i?(r')}  =E(a)|gyl(^)|E(0)|?B(7■,)}  (46) 
and 

E (t){&4 (»*)&? (»•')}  =  E(6)j<3U(r)}  E(t)|9B(r')}  (47) 
respectively.  From  Eq.  (37)  and 

E(6)  { Qb (r') }  =  *52  Cj  Pd  (r'>  Ri)>  (48) 

i=i 

we  can  obtain  the  following  result, 


E(e){gA(r)9B(r')} 

\n—j 

=  R4C  f  dRp$(r,R)pa(r',R)  qin(R).  (51) 

J  OO 

Next,  consider  the  case  of  n  ^  j  for  step  (c).  Prom 
Eq.  (49),  we  have 


E(c){<m(*’)4b(0} 


n^j 


N  N 


=  (EE^OPd(r,Hn)pf(r-,,f?i) 


n= 1 j= 1 
N  N 


(<=) 


n^j 


Cj}Pd(r’Rn)Pd(r',Rj) 

n=l j— 1 

= ic  E  Ep^r’ Rn ) p*  (r'’ 

n= 1 j— 1 


(52) 


n-=Aj 


Again,  similarly  as  in  Eq.  (39),  step  (d)  applied  to 
Eq.  (52)  now  gives 

E(d){^(r)Mr,)}|  \  dR  f  dR' UA\r,R) 

^  J  I  n^j  doo  Joo  t 

xriB,(-''.H')i;Eprft.,ftj(«.R'lW)}|  •  <53> 

n— 1 j—l  J  I n^j 


N  N 

E  (6)  {^4  (O'Mr')}  =  EE^OPd(r,Hn)pf(r',fii). 

n— 1 j= 1 

(49) 

In  order  to  average  over  {|n}  and  {£/}  for  fixed 
qin(r)  in  step  (c),  we  must  consider  two  cases,  de¬ 
noted  by  E(c){ofA(^)^B(^,)}|n=j  when  n  —  j  and,  by 
E (c){QA(r)qB(rf)}\n^j  when  n  ^  j.  For  the  double 
sum  over  n  and  j  in  Eq.  (49),  there  are  N  terms  with 
n~  j.  Averaging  these  N  terms  in  Eq.  (49),  yields 


E(o) 


n~  3 


Based  on  the  work  by  Barrett  et  al  [19],  it  can  be 
shown  that 


N  N 


E{qin(R)qin(R')\N }  =  £  E  WknMR’  R'\ N) 

k  n^j  n=i j=l 

(54) 


Substituting  Eq.  (54)  into  Eq.  (53),  we  have 

E(d) { qA (r)qB (r') } I  =Zt[  dR  [  dR' l Pd(r, R) 

{  *  | n-£j  Joo  Joo  L 

xpB(r',R')E{qin(R)qin(R')\N}\  }.  (55) 
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Now  it  is  easily  shown  that  step  (e)  applied  to  Eq.  (55), 
to  obtain  the  average  over  iV,  gives 


E(*){«i4(r)&(p')} 

\n^j 

=  tt[ W dR'pfrr,  R)Pj(r\ R')E{qin(R)qin(R')} 

J  OO  J  OO 

=  tc[  dR[  dR'pA(r,R)p*(r',R')Rin(R,R')\  .  (56) 

J  oo  J  oo  ln*» 


To  replace  R,n(.R,  in  Eq.  (56)  with 

R in(R,R'),  we  invoke  Eqs.  (26)-(27)  again.  Thus 
Eq.  (56)  becomes 
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E(e){gA(r)gB(p')} 

'  n&j 

=  U  f  dR  f  dR'p^(r,R)p^(r',R')Rin(R,R') 

J  oo  «/oo 

-£C  f  dRp$  (r,  R)pj  (r1,  R)  qin(R).  (57) 

J  OO 

Adding  Eqs.  (51)  and  (57),  the  cross  correlation  of 
qA(r)  and  qB(r)  is  given  by 


FIG.  3.  Two  point  processes  following  multiple  ampli¬ 
fied  point  processes  subsequent  to  selection  of  a  random 
point  process. 

amplification  stages  in  each  of  the  two  paths  A  and  B, 
as  illustrated  in  Fig.  3.  In  doing  so,  it  is  important  to 
note  that  there  are  no  random  variables  in  Eq.  (59), 
the  cross  covariance  is  expressed  only  in  terms  of  mean 
values  and  the  cross  covariance  of  |n  and  £n. 


Rab(p,0  =  k€C  [  dRp${r,R)p$(r',R)qin(R) 

J  OO 

+fC  /  dR  j  dR'p2(r,R)p$(r',R')Rin(R,R'),  (58) 

J  oo  J  oo 

which  is  a  general  expression  for  the  non-stationary 
cross  correlation  of  two  amplified  point  processes  fol¬ 
lowing  a  stochastic  input  selection.  The  cross  covari¬ 
ance  of  qA{r)  and  qB(r)  is  given  by  subtracting  the 
product  of  means  from  the  cross  correlation.  From 
Eqs.  (42),  (44)  and  (58),  finally,  we  have 


1.  Transfer  Function  of  Mean 


We  approach  the  generalization  of  Eq.  (59)  by  first 
considering  the  general  expression  Eq.  (43)  for  the 
mean  of  an  amplified  point  process  with  any  input 
point  process.  Without  loss  of  generality,  we  define  a 
transfer  function  iJ(-),  an  integral  operator,  in  terms 
of  Eq.  (43),  as 


KAB(r,r')=Kw  /  dRpj(r,R)pa(r',R)qin(R) 

J  OO 

+U  [  dR  [  dR'p^(r,R)p^(r',R')Kin(R,R'),  (59) 
J  OO  J  OO 

which  is  the  desired  result.  For  the  case  illustrated 
in  Fig.  2,  the  expressions  given  by  Eqs.  (58)  and  (59) 
show  that:  a)  the  correlation  in  qA(r)  and  qB(r )  is 
proportional  to  the  cross  correlation  of  the  binomial 
random  variables  £n  and  (n  describing  the  point  se¬ 
lection  as  given  by  the  first  term  on  the  right-hand 
side  of  Eqs.  (58)  and  (59);  and,  b)  any  correlation  in 
the  random  source  qin{r)  is  transferred  to  the  outputs 
through  paths  A  and  B  as  shown  by  the  second  term. 

D.  Cross  Covariance  Following  Multiple  Amplified 
Point  Processes 

We  now  generalize  the  result  of  Eq.  (59)  derived  in 
the  above  section  to  an  arbitrary  number  of  cascaded 


H{qin(R))=  j  dRpd(r,R)qin(R),  (60) 

which  describes  the  propagation  of  the  mean  of  a 
point  process  passing  through  an  amplification  stage. 
As  shown  in  Fig.  4(a),  given  an  amplification  stage 
with  transfer  function  H(-) ,  the  mean  of  the  ampli¬ 
fied  point  process  qout{r )  resulting  from  any  point 
process  input  qini^)  can  be  uniquely  determined,  i.e., 
qout(r)  =  H(qin(R)),  where  R  and  r  represent  the  po¬ 
sition  vectors  of  input  and  output  quanta,  respectively. 
The  mathematical  operation  in  Eq.  (60)  that  maps  the 
mean  of  input  at  R  to  the  output  at  r  is  a  convolution 
of  the  input  qin(R)  with  the  function  pd(r,R).  For 
convenience,  denote  to  represent  the  “shift  variant” 
convolution  integral  such  that 

H^qin(R) )  =  J  dRpd{r,R)qin(R) 

=  Pd{r,R)  *v  Qin(R),  (61) 
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qin(R) 


H(  •  ) 

(a) 


Qout  (v ) 


Qini^fy 


*v  Pd(r,R ) 

(b) 


Qout{^) 


(c) 

FIG.  4.  Transfer  functions  of  mean  for  an  amplified 
point  process. 


where  the  subscript  v  of  the  symbol  *  is  used  to 
make  a  distinction  between  shift  invariant  and  vari¬ 
ant  systems.  Fig.  4(b)  is  the  block-diagram  repre¬ 
sentation  of  Eq.  (61).  From  the  relationship  between 
function  pd{r ,  R)  and  the  probability  density  function 
pr^(r  —  H|.R)  given  by  Eq.  (31),  we  have 

H(qin(R)J=  J  dR pr ±(r  -  R\R)kn(R)qin(R) 

=  prA(r  -  R\R)  *v  ( kn(R )  qin(R)) ,  (62) 

which  is  an  equivalent  expression  of  the  mean  for  an 
amplified  point  process.  As  illustrated  in  Fig.  4(c),  an 
input  point  process  qin(R )  is  first  amplified  with  mean 
gain  kn{R)  and  then  scattered  with  the  density  func¬ 
tion  pr^(r  -  R\R).  This  result  can  be  found  in  [13] 
developed  by  Rabbani  et  al  from  the  view  of  multi¬ 
variate  moment-generating  functions. 

If  the  amplified  point  process  is  shift-invariant,  then 
the  mean  gain  and  density  function  are  independent 
of  position  R ,  where  kn(R)  =  k  and  pr^(r  -  R\R)  — 
pr^(r  -  R).  In  this  case,  Eq.  (62)  becomes 

H(qin(R)Sj=  J  <LRprA(r  -  R)kqin(R) 

=  kpr^(r)*qin(r).  (63) 


2.  Cross  Covariance 

Now  consider  the  case  illustrated  in  Fig.  3,  where 
there  are  L  and  M  cascaded  amplification  stages  in 
paths  A  and  B,  respectively.  The  output  of  each  am¬ 
plification  stage  forms  a  virtual  input  to  the  next. 


First,  using  an  approach  similar  to  the  procedure 
discussed  previously,  we  calculate  the  cross  covariance 
of  the  outputs  for  the  case  of  two  cascaded  amplifi¬ 
cation  stages  in  each  of  the  paths  A  and  B.  In  other 
words,  we  may  make  the  calculation,  for  that  case,  by 
considering  the  five  steps  shown  in  the  previous  sec¬ 
tion.  Then,  by  the  method  of  induction,  the  cross  co- 
variance  for  the  output  of  the  system  shown  in  Fig.  3 
can  be  obtained.  However,  this  procedure  becomes 
tedious  due  to  the  conditional  average  over  multiple 
amplification  processes.  Rather,  cascaded  amplifica¬ 
tion  stages  in  each  path  can  be  thought  of  as  an  entire 
amplified  point  process  if  we  consider  only  linear  sys¬ 
tems  where  the  spatial  pattern  of  each  amplification 
stage  depends  only  on  the  position  of  the  primary  in¬ 
teraction.  To  obtain  the  statistical  characteristic  for 
the  entire  amplified  point  process,  consider  the  calcu¬ 
lation  of  the  mean  of  the  output  point  process  qA(r) 
by  cascading  Eq.  (61)  with  the  L  amplification  stages 
in  path  A.  In  Fig.  5,  the  mean  qA(r)  of  the  output  at 
r  is  the  convolution  of  qAh_  1(RAl~1)  with  the  func¬ 
tion  pAL(r:RAL~1)  and  then  of  qAL_2{RAL~2)  with 
the  function  pAh~x  {RAl-\RAl~2),  •*,  finally  of  in¬ 
put  at  R  with  the  function  pAl(RAl,  R).  This  lead 
to 

QA(r)  =p$L(rtRAL-')  *vqAL-i(RAL-') 

=  pA*{r,  RA^)  *„  [PA-{RAL-\  RAl~2)  *vqAL_,(RAL-2)} 

=  PAl(t,RAl-')  *v---*vPAi(rAi,R)  *v  Qa0{R) 

—Pd(riR)  *v  Qa0{R),  (64) 

where  we  define  pA(r,R)  as  the  statistical  character¬ 
istic  function  for  the  cascade  of  all  stages  along  path 
A  given  by 

pA(r,R)^pAl-(r,RA^)  *„  ■  •  •  *„  Pa'(ra\r). 

(65) 

Now  we  can  say  that  the  mean  of  the  output  q(r)  is 
the  convolution  of  the  mean  of  the  input  at  R  with 
the  function  pA(r ,  R ). 

Similarly,  for  path  B,  we  obtain  that  the  mean  of 
qB(r)  is  given  by 

qB(r)=p$(r,R)*vqBo(R),  (66) 

where  point  process  qe0  (r)  is  the  output  of  the  point- 
selection  process  for  path  B  (see  Fig.  3)  and 

pf  (r,  R)^M(r,  RBm-1)  *„•••*„  PBd\RBl ,  R)  ■ 

(67) 
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N 


Qaq(R) 


*v  Pd(r  ’  R) 


U(r) 


FIG.  5.  Transfer  function  of  mean  for  multiple  amplified  point  processes. 


This  result  shows  that  multiple  cascaded  amplified 
point  processes  can  be  expressed  as  a  single  amplified 
point  process  with  the  appropriate  statistical  charac¬ 
teristic  function  for  the  purpose  of  describing  transfer 
of  mean  value.  Since  calculation  of  the  cross  covariance 
in  Eq.  (59)  involves  integrals  in  which  amplification  is 
represented  in  terms  of  the  statistical  characteristic 
only,  we  conclude  that  the  cross  covariance  for  qA  (r) 
and  qB  (r)  is  also  given  by  Eq.  (59),  where  pj(r,R) 
is  described  by  Eq.  (65)  and  pB  (r,  It)  is  described  by 
Eq.  (67). 


E.  Important  Special  Cases 
1.  Doubly  Stochastic  Poisson  Impulses  Source 


where  R^  =  K#  +  £  C  is  the  cross  correlation  of  ran¬ 
dom  binomial  variables  £n  and  Similarly,  substi¬ 
tuting  Eqs.  (68)  and  (69)  into  Eq.  (59),  after  some 
algebraic  manipulations,  the  cross  covariance  becomes 

KAB(r,r')  =  R?<  [  dRpj{r,R)p%{r',R)b{R) 

J  OO 

+H  [  dR  f  dR’p$(r,R)p$(r',R')Kb(R,R'),  (71) 

which  is  the  result  of  Fig.  2  for  the  input  process  of 
doubly  stochastic  Poisson  impulses.  Eqs  (70)  and  (71) 
are  given  to  show  the  relationship  to  Barrett’s  work, 
but  are  not  required  to  obtain  the  following  special 
cases. 


In  medical  imaging  systems,  the  incident  x  rays  at 
the  detector  form  a  quantum  image  which  is  a  dis¬ 
tribution  of  quanta,  and  can  be  expressed  as  a  sample 
function  of  a  spatial  point  process  qin{r )  as  in  Eq.  (15). 
Barrett  [19]  has  shown  that  such  a  point  process  may 
be  described  in  terms  of  doubly  stochastic  Poisson  im¬ 
pulses  with  intensity  process  6(r),  which  he  calls  the 
random  input  fluence.  In  x  ray  imaging,  for  instance, 
random  changes  of  the  radiation  source  makes  b(r) 
spatially  random. 

For  the  input  point  process  qin(r)  of  doubly  stochas¬ 
tic  Poisson  impulses,  Barrett  shows  that  the  mean  is 
given  by 

qin(r)  =  b(r),  (68) 

where  b(r)  is  the  mean  of  random  input  fluence,  and 
the  autocovariance  is 

Km(r,0  =  b(r)S(r  -  r')  +  Kb{r,r'),  (69) 

where  Kb(r,r')  is  the  autocovariance  of  b(r).  Sub¬ 
stituting  Eqs.  (68)-(69)  into  Eq.  (29),  we  obtain  (see 
Fig.  1) 

K AB(r,r')  =  R^  b(r)6(r  -  r')  +  £(Kb(r,r'),  (70) 


2.  Shift-Invariant  System  with  Multiple  Amplified  Point 
Processes 

If  the  system  shown  in  Fig.  3  is  shift-invariant,  which 
requires  that  the  mean  gain  and  the  probability  den¬ 
sity  function  for  each  amplification  stage  be  indepen¬ 
dent  of  position,  then  the  propagation  of  the  mean  of 
the  input  point  process  is  shown  in  Fig.  6.  From  Fig.  6 
and  Eq.  (63),  it  is  easily  shown  that 

qA(r)  =  kA  pr^(r)  *  qAo{r ) 
qB(r)  =  kB  pr£(r)  *  qBo{r ) 

where  the  probability  density  functions  pr^(r)  and 
pr?  (r)  for  the  entire  amplified  point  processes  along 
paths  A  and  B,  respectively,  are  expressed  as  convolu¬ 
tions  of  the  density  functions  of  sub-stages  in  its  path, 
i.e., 

pr-4  (r)  =  p4l  (r)  *  pr^”1  (r) *  •  •  •  *  pr£  (r) 
prf  (r)  =  pr|M  (r)  *  pr|M~’  (r)  *  •  •  •  *  pr|J  (r) 

(73) 
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FIG.  6.  Shift-in variant  system  with  multiple  amplified 
point  processes. 

and  the  mean  gains  kA  and  kB  for  the  entire  ampli¬ 
fied  point  processes  of  paths  A  and  B,  respectively,  are 
given  by 


< 


kA = n 

i= 1 
M 

kB = n  *Bi 

3= 1 


(74) 


From  Eq.  (31),  we  obtain  the  function  p^(r,  R)  for  the 
entire  amplified  point  process  of  path  A,  i.e., 

pA(r,R)  =  kApvl(r-R).  (75) 


Similarly,  we  have 


covariance  K ^(r,/).  For  the  wide-sense  stationary 
conditions,  the  input  process  must  be  stationary  in 
the  wide  sense  and,  the  amplification  processes  must 
be  shift-invariant  with  uniform  mean  gains  of  quanta 
in  an  infinite  imaging  plane.  Thus,  under  WSS  condi¬ 
tions,  let  qin(R)  =  qin  be  constant,  and  Kin(R,R')  = 
Kin(R-R').  Moreover,  the  functions  pA(r,R)  and 
p%(r',R')  in  Eq.  (59)  are  replaced  with  the  products 
of  their  corresponding  constant  mean  gains  and  shift- 
invariant  density  functions  as  in  Eqs.  (75)  and  (76), 
respectively.  We  obtain 

K AB{r,r')  -  K AB(r  -  r#) 

=  Ki(kAkB  qin  ( dHpr4(r-jR)pr|(r'-i?) 

J  OO 

+ZCkAkB[  <mj  dR'prjjr  -  R)pr|(r'  -  R')Kin(R- R'), 

J  OO  J  OO 

(78) 


By  changing  the  integral  variable  R  such  that  r'  — 
r'  -  R ,  the  integral  of  the  first  term  in  Eq.  (78)  can 
be  rewritten  as 

[  dr'  pr£  (r-r1  +  r')  prf (r') 

J  OO 

which  is  the  correlation  integral  of  two  functions 
pr^(r)  and  pr£(r)  over  r',  denoted  in  short  form 
by  pr^(r)  *pr£(r),  where  r  =  r  -r'.  Moreover,  for 
the  double  integral  of  the  second  term  in  Eq.  (78),  we 
perform  the  changes  of  t"  =  r'  —  R!  in  the  R!  integral 
and  r"'  +  t"  =  r  -  R  in  the  R  integral,  yielding 


Pd  (r,  JR)  =  kB  pr|(r  -  R).  (76) 

Substituting  Eqs.  (75)  and  (76)  into  Eq.  (59),  the  cross 
covariance  function  of  gUO")  and  <7#(r)  is  given  by 

K  AB(ry) 

=  K ^kAkBf  dR pr£(r - R)  pr|(r' -  R)  qin (R) 

J  OO 

-K<  kAkB[  dRf  dR'pri(r-R)vrl(r'-R')Kin(R,R') 

J  CO  J  OO 

(77) 


[  dr'"  f  dr"pr^(r'"  +  r")  prf(r")  Kin(r-r'-r'") 

J  OO  L  ^  OO 


The  integral  over  r"  is  the  correlation  integral  of 
prd(r'")  and  pr£(r"'),  and  the  integral  over  r'"  is 
the  convolution  integral  of  the  correlation  integral  and 
K in  (t)  .  At  this  point  we  simplify  our  notation  and  let 
pr (t)  =  pr^(r).  Therefore,  Eq.  (78)  becomes 


K Ab(t)  =  K^kAkBqin  prA(r)  *prB(r) 

+£ (  kA kB  jpr"4 (t)  *  prB (r ) j  *  Kin (r ) ,  (79) 


which  is  the  desired  result  for  the  case  when  the  system 
is  shift  invariant. 

3.  Cross  Covariance  and  Cross  Spectral  Density  Under 
WSS  Conditions 

If  the  outputs  of  point  processes  £4(1*)  and  Qb{t) 
in  Fig.  2  are  wide-sense  stationary  (WSS),  we  can  de¬ 
scribe  the  correlation  between  two  paths  A  and  B  in 
the  frequency  domain  by  their  cross  spectral  density 
[3]  which  is  equal  to  the  Fourier  transform  of  the  cross 


which  is  the  desired  expression  of  the  cross  covariance 
for  wide-sense  stationary  gU(^)  and  £s(r). 

The  cross  spectral  density  for  paths  A  and  B  is  de¬ 
fined  as  the  Fourier  transform  of  the  cross  covariance 
Kab(t),  given  by 

NPS^bM  =  ^{k^b(t)} 

=  K? c  F1  P  <?*„  .T7  {  pr 4  (t)  *  prB  (t)  } 

+^C^fcB^|[prAW*PrB(r)]  *Kin(r)}.  (80) 
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FIG.  7.  Schematic  of  the  parallel  cascade  model  used  to  determine  the  reabsorption  of  characteristic  K  x  rays  in  a 
radiographic  screen. 


Then  the  final  result  for  the  cross  spectral  density  un¬ 
der  WSS  conditions  is 

NPS  AB{v)  =  kAkB  TV)TB»  [£CNPSin(i/)  +  K^qin 

(81) 

where  NPSf n(u)  is  the  NPS  of  the  input  point  process 
(the  Fourier  transform  of  Kjn(r)),  TA(ts)  and  T B(u) 
are  Fourier  transforms  of  prA(r)  and  prB(r)  respec¬ 
tively.  When  paths  A  and  B  represent  a  cascade  of 
multiple  amplified  point  processes,  the  mean  gains  kA 
and  kB  are  the  product  of  mean  gains  along  each  path 
as  given  by  Eq.  (74),  and  pr^r)  and  prB(r)  rep¬ 
resent  the  cascaded  probability  density  functions  for 
each  path  as  given  by  Eq.  (73).  Taking  the  Fourier 
transform  of  both  sides  of  Eq.  (73)  shows  that  TA(is) 
and  Tb  (is)  represent  the  product  of  all  scatter  transfer 
functions  along  each  path,  as  given  by 

T»-nT*M 

I1  ■  w 

tb{v ) = nTB,» 

j= i 

where  TAi  (is)  is  the  Fourier  transform  of  the  probabil¬ 
ity  density  function  of  the  zth  sub-stage  in  path  A.  The 
probability  density  functions  always  have  unity  area, 
and  hence  the  transfer  functions  T A(v)  and  TB(is) 
always  have  a  value  of  unity  at  v  —  0.  While  the 
transfer  functions  are  complex  in  general,  as  shown 
in  Appendix  B,  the  sum  of  any  cross  term  pair, 
NPS^b(^)  +  NPS ba(v),  will  always  be  real  only. 


III.  APPLICATIONS 

Equation  (81)  is  the  general  expression  for  the  WSS 
cross  spectral  density  of  two  parallel  cascaded  ampli¬ 
fied  point  processes  descending  from  a  single  input 
point  process.  It  is  in  a  particularly  simple  and  con¬ 
venient  form  for  application  in  a  linear-systems  model 
as  it  is  expressed  solely  in  terms  of  the  NPS  and  mean 
value  of  the  input  point  process,  the  selection  proba¬ 
bilities  and  corresponding  cross  covariance,  and  mean 
gains  and  mean  scatter  transfer  functions  for  each 
path.  In  this  section,  this  result  is  used  in  a  description 
of  characteristic  reabsorption  in  a  radiographic  screen. 


A.  Application  to  Reabsorption  of  Characteristic 
X  Rays  in  a  Radiographic  Screen 

We  examine  here  the  effects  on  image  noise  of  fluo¬ 
rescence  reabsorption  in  a  radiographic  screen,  where 
light  is  emitted  at  both  the  primary  photo-electric  in¬ 
teraction  site  and  at  the  reabsorption  site.  This  prob¬ 
lem  was  studied  previously  by  Metz  and  Vyborny  [34] 
using  a  relatively  sophisticated  statistical  analysis.  We 
show  that  the  same  result  can  be  obtained  using  a  sim¬ 
pler  linear  transfer-theory  model  that  includes  parallel 
cascades  and  the  cross  spectral  density  derived  in  the 
previous  section. 

Figure  7  illustrates  a  “flow  diagram”  showing  the 
sequence  of  events  leading  to  light  production  in  the 
Metz- Vyborny  model.  WSS  conditions  are  assumed 
throughout  so  that  each  position  in  this  diagram  rep¬ 
resents  an  intermediate  step  between  input  and  output 
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characterized  in  terms  of  a  two-dimensional  distribu¬ 
tion  of  quanta  (points)  q(r).  The  processes  included  in 
Fig.  7  are  based  on  three  “elementary  processes”  (see 
Appendix  A)  in  the  serial  cascades  plus  branch-points 
that  give  rise  to  the  parallel  cascades. 

Several  simplifying  assumptions  are  made  in  order 
to  be  consistent  with  Metz  and  Vyborny  [34]  and  with 
earlier  work  of  Rossmann  [35,36].  They  include:  a)  in¬ 
cident  x  rays  are  assumed  monoenergetic;  b)  differ¬ 
ences  in  light  emission  due  to  different  x  ray  interaction 
depths  are  ignored;  and  c)  only  photo-electric  interac¬ 
tions  are  considered.  Metz  and  Vyborny  also  ignored 
the  statistical  nature  of  light  generation  in  the  screen. 

At  the  input  of  the  model  in  Fig.  7,  a  uniform  x-ray 
distribution  consisting  of  q  quanta/mm2,  each  with  en¬ 
ergy  Ex,  is  incident  on  the  radiographic  screen.  These 
quanta  are  Poisson  distributed,  and  hence  have  an  as¬ 
sociated  NPS  given  by  NPS(i')  =  q  [37].  A  fraction  rj 
of  these  incident  quanta  will  result  in  a  photo-electric 
interaction  in  the  screen.  Selection  of  these  events 
is  represented  as  a  stochastic  selection  or  binary  gain 
stage,  where  gain  is  represented  by  a  random  variable 
rj  that  can  have  a  value  of  0  or  1  only  and  mean  of  rj. 
The  output  from  this  gain  stage  is  the  two-dimensional 
distribution  of  photo-electric  events  in  the  screen. 

As  described  by  Metz  and  Vyborny,  there  are  three 
possible  sequences  of  events  whereby  light  can  be  gen¬ 
erated  for  each  photo-electric  interaction:  1)  absorp¬ 
tion  of  the  primary  x-ray  photon  at  the  primary  in¬ 
teraction  site  without  emission  of  a  characteristic  K  x 
ray;  2)  absorption  of  the  primary  x  ray  accompanied 
by  emission  of  a  K  x  ray;  and  3)  reabsorption  of  the 
K  x  ray  at  a  remote  location.  These  three  sequences 
correspond  to  paths  A,  B  and  C  in  Fig.  7. 

Path  A  describes  the  emission  of  light  at  the  primary 
interaction  location  when  no  K  x  ray  is  produced.  For 
each  photo-electric  interaction,  there  is  a  probability 
qu  that  a  K  x  ray  will  be  generated,  and  therefore  a 
probability  (1  -  qoS)  that  a  K  x  ray  is  not  generated 
where  q  is  the  probability  that,  when  an  incident  pho¬ 
ton  interacts  in  the  screen,  it  undergoes  a  K-shell  inter¬ 
action,  and  c 0  is  the  fluorescent  yield  of  K-shell  photo¬ 
electric  interactions.  This  branching  is  represented  in 
Fig.  7  as  the  diamond-shaped  “Bernoulli  branch”.  It 
is  to  be  interpreted  as  a  Bernoulli  trial  [2,37]  that, 
for  each  interaction,  determines  the  outcome  “yes”  or 
“no”  where  “yes”  is  obtained  randomly  with  a  prob¬ 
ability  qu,  and  “no”  otherwise.  If  a  K  x  ray  is  not 
produced,  corresponding  to  path  A  in  Fig.  7,  it  is  as¬ 
sumed  that  the  incident  x-ray  energy  Ex  is  absorbed 
locally  producing  the  number  m  optical  quanta  (the 
gain  factor  m  is  assumed  to  be  proportional  to  the  ab¬ 
sorbed  energy)  which  will  be  emitted  from  the  screen. 
Metz  and  Vyborny  ignore  the  statistical  nature  of  light 
emission,  and  hence  m  is  modeled  as  a  deterministic 
gain  factor.  Note  that  the  gain  factor  m  only  describes 
generation  of  the  number  of  light  quanta  that  are  emit¬ 
ted  from  the  screen. 


Path  B  describes  light  emission  at  the  site  of  the 
photo-electric  interaction  when  a  K  x  ray  is  emitted 
(which  may  or  may  not  be  reabsorbed).  In  this  case, 
the  energy  Ek  —  ^Ex  is  carried  away  in  the  K  x  ray, 
and  the  remaining  energy  Ex  —  Ek  is  deposited  at  the 
primary  interaction  site.  Thus,  only  (1  —  -0)m  optical 
quanta  are  emitted  at  the  primary  interaction  site  for 
each  photo-electric  interaction  where  ^  =  1 Ek/Ex  and 
Ek  &  59.3  kev  for  tungsten  in  the  calcium  tungstate 
screens. 

Path  C  describes  the  light  emitted  from  the  screen  at 
a  remote  site  due  to  reabsorption  of  the  K  x  ray,  where 
fK  is  the  probability  of  reabsorption  somewhere  in  the 
screen  for  each  photo-electric  interaction  producing  a 
K  x  ray.  The  location  of  reabsorption  is  random,  but 
the  point-spread  function  Pk{v)>  which  has  unity  area, 
describes  the  probability  density  that  the  K  x  ray  is 
reabsorbed  at  a  distance  r  from  the  photo-electric  in¬ 
teraction  site.  The  process  representing  this  random 
relocation  of  the  K  x  ray  is  a  quantum  scatter  stage 
as  described  in  Appendix  A.  At  the  reabsorption  site, 
the  K  x  ray  is  converted  to  optical  quanta  with  a  con¬ 
version  factor  tym.  Events  are  selected  for  both  paths 
B  and  C  for  every  “no”  event  in  the  Bernoulli  branch. 
We  call  the  point  of  separation  of  paths  B  and  C  a 
“cascade  fork.” 

Due  to  geometrical  spread  and  possibly  light  scatter 
in  the  screen,  optical  quanta  are  distributed  spatially 
with  a  point  spread  function  (normalized  to  unity  area) 
given  by  p0(r ),  Metz  and  Vyborny  ignore  the  statis¬ 
tical  nature  of  light  scatter  and  this  redistribution  of 
light  is  represented  as  a  linear  filter  (convolution)  with 
a  kernel  p0(r). 

The  total  light  emitted  from  the  screen  is  therefore 
the  sum  of  contributions  from  each  path,  resulting  in 

q0  =  \To(0)\[qA  +  qB+  qc\  (83) 

quanta  per  unit  area,  where  we  have  used  Eq.  (A9) 
for  the  output  linear  filter  in  Fig.  7,  and  the  function 
T0{u)  is  the  Fourier  transform  of  p0(r).  The  contri¬ 
butions  from  each  path  can  be  obtained  by  cascadeing 
the  elementary  processes  (see  Appendix  A)  included 
in  each  path  and  considering  the  outcome  probability 
of  the  Bernoulli  branch.  Then,  we  have 

qA  -  m(l  —  qco)r]  <?,  (84) 

qs  =  (1  -  ^mqcupq,  (85) 

qc  =  ipmfKWn  q •  (86) 

Substituting  Eqs.  (84)-(86)  into  Eq.  (83),  we  can  ob¬ 
tain 

qQ  =  qrjm[l  -  qunp(  1  -  /*:)]|To(0)|.  (87) 

The  NPS  of  light  emitted  from  the  screen,  denoted 
by  NPS0(u),  is  therefore  given  by  (see  Appendix  B) 
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t  APP:  Amplified  Point  Process 
FIG.  8.  The  Bernoulli  branch  with  amplified  point  pro- 


^APP:  Amplified  Point  Process 
FIG.  9.  The  cascade  fork  with  amplified  point  processes. 


NPSo(*/)  =  [nPSaH  +  NPSs(i')  +  NPSc(^) 
+NPSabM  +  NPS  BA(v)  +  UPSacM 
+NPSca(^)  +  NPSbcM  +  NPScbM]  \T0(v)\2  (88) 

consisting  of  the  NPS  from  each  of  the  paths  A,  B 
and  C  plus  corresponding  cross  terms  as  described  be¬ 
low  for  the  parallel  paths  with  “Bernoulli  branch”  and 
“cascade  fork”  selection  processes. 


1.  Bernoulli  Branch 

The  Bernoulli  branch  with  amplified  point  processes 
is  illustrated  in  Fig.  8.  Each  quantum  in  the  in¬ 
put  point  process  is  selected  for  path  A,  denoted  by 
“yes”,  when  |n  =  1  and  for  path  B,  denoted  by  “no” 
when  £n  =  0.  The  Bernoulii  branch  is  a  special 
case  of  the  point  selection  process  described  in  Fig.  2 
where  the  two  binomial  random  variables  are  related 
to  C  =  (1  -  £).  This  results  in  the  cross  covariance  of 
these  random  variables  given  by 

K"  =E{|„Cn}  -E{|n}E{c„}  =  -£<  =  -£(1- A 

(89) 

The  cross  spectral  density  following  amplification  is 
therefore  given  by  Eq.  (81)  as 

NPS  AB(v) 

=  |(1  -  t)kAkB  tV)Tb»  [NPSin(*0  -  gin]  (90) 

showing  that  there  is  correlation  between  paths  A  and 
B  only  if  quanta  in  the  input  image  are  statistically 
correlated.  That  is,  when  NPSf„(i/)  —  qin  #  0.  If  the 
quanta  are  uncorrelated  and  NPSin(k')  =  qin>  there  is 
no  cross  term. 


2.  Cascade  Fork 

The  cascade  fork  with  amplification  is  shown  in 
Fig.  9  where  every  quantum  in  the  input  is  selected 
for  both  paths  A  and  B.  This  again  is  a  special  case  of 
the  general  point  selection  process  described  in  Fig.  2 
where  £n  =  (n  =  1,  modeled  as  deterministic  unit  fac¬ 
tors.  The  cross  covariance  of  £n  and  £n,  is  therefore 

K«  =  E{fnCn}  -  e{|„}e{c„}  =  0,  (91) 

and  the  cross  spectral  density  term  for  the  cascade  fork 
based  on  Eq.  (81)  is  given  by 

NPS ab{v)  =  kAkB  TV)TB>)NPS»»,  (92) 

which  is  always  non-zero  if  the  input  is  a  random  point 
process,  and  therefore  there  is  always  a  cross  term  be¬ 
tween  paths  A  and  B. 

B.  Degradation  of  the  NPS  due  to  Reabsorption 

The  NPS  in  the  distribution  of  optical  quanta  from 
each  path  in  isolation  is  obtained  by  cascading  appro¬ 
priate  combinations  of  the  elementary  processes  de¬ 
scribed  in  Appendix  A,  giving 


NPSaM  =m2(l  -<rw)7?g, 

(93) 

NPSbM  =  (1  -  ip)2m2<;uT]q, 

(94) 

NPScM  =  4>2m?fk<;uriq, 

(95) 

where  we  have  used  NPS(i^)  =  q  since  the  input  quanta 
are  Poisson  distributed.  Quanta  in  the  input  image  are 
statistically  uncorrelated,  therefore  there  is  no  cross 
term  between  paths  separated  by  the  Bernoulli  branch. 
That  is, 

NPSj4B(i/)-NPSbaW  =  0, 

WPSacW  =  NPS  caW)  =  0- 


Appendix  I:  Parallel  Cascades 


24 


FIG.  10.  Illustration  of  Tfe(i/),  the  MTF  corresponding 
to  reabsorption  in  the  film-screen  system  (adapted  from 
Metz  and  Vyborny). 
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FIG.  11.  Illustration  of  F(i/)  from  Eq.  (100),  showing 
contributions  from  the  uncorrelated  and  correlated  com¬ 
ponents  (adapted  from  Metz  and  Vyborny). 
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The  cross  terms  NPS bc{v)  and  NPS cbM,  based 
on  Eq.  (92),  are  given  by 

NPS Bc(y)  =  (1  -  '<P)'4’m2fKT*K{y)<;w'q q,  (96) 

and 

NPScb(^)  =  (1  -  ip)ipm2fKTidt')‘;vvqi  (97) 

respectively,  where  Tk{v)  is  the  characteristic  transfer 
function  describing  the  reabsorption  probability  den¬ 
sity  in  terms  of  spatial  frequencies  and  is  equal  to  the 
Fourier  transform  of  the  reabsorption  PSF,  p/c(r).  It 
is  known  that  the  sum  of  two  complex  conjugates  is 
equal  to  two  times  the  real  part.  This  leads  to 

NPSsc(i')  +  NPScb(^) 

=  2qri<;u)(l  -  •0)i/>m2/ArRe(Tie(i/)),  (98) 

which  is  always  real,  where  Re{  }  denotes  the  real  part 
of  a  complex  quantity. 

Combining  the  above  results  gives  the  NPS  for  the 
output  optical  image  quanta,  including  the  effect  of 
the  redistribution  of  light  in  the  screen,  as 

NPSo(i')  =  qr/m 2  [(1  -  qu)  +  ?w(l  -  ip)2  +  qcjfxip2 

+2w/kV’(1  -  V')Re{Tif(^)}]|To(^)|2  (99) 

which  is  the  Metz-Vyborny  result  for  the  NPS  of  light 
emitted  from  the  screen. 

Metz  and  Vyborny  used  Eq.  (99)  to  describe  the  ef¬ 
fect  of  reabsorption  in  a  Dupont  Par  Speed  calcium 
tungstate  screen.  Assuming  a  constant  film  density 
(fixed  light  output  q0),  the  effect  of  reabsorption  is 
obtained  by  considering  the  NPS  just  above  and  be¬ 
low  the  K-edge  of  tungsten  at  approximately  68.5  keV . 
Above  the  K-edge  reabsorption  takes  place  and  q0  is 
given  by  Eq.  (87).  Below  the  K-edge,  both  the  light 
output  and  the  NPS  are  determined  by  setting  ip  =  0 
and 


=  0.  Therefore,  the  ratio  of  the  NPS  just  above  the 
K-edge  to  just  below,  normalized  to  fixed  total  light 
output,  is  given  by  T(i/)  where 

I»  = 

(l  -  w)-Kw(l -ip)2+qufKip2+2<;ufK,4>Q-  -^)Re  (Tk^)) 

1  -  quiip(l  -  fa) 

(100) 

Figure  10  illustrates  T k(v)  as  used  by  Metz  and 
Vyborny.  The  corresponding  degradation  in  the  NPS, 
r(u),  is  shown  in  Fig.  11  obtained  using  values  listed  in 
Table  1.  The  value  of  NPS(n-)  is  decreased  by  approx¬ 
imately  10%  at  low  frequencies,  increasing  to  within  a 
few  percent  of  the  uncorrelated  value  for  frequencies 
above  approximately  2  cycles/mm.  While  these  re¬ 
sults  are  specific  to  a  calcium  tungstate  screen  which 
has  limited  use  at  present,  corresponding  results  for 
newer  screens  and  other  imaging  systems  can  be  ob¬ 
tained  using  the  same  formalism.  The  transfer-theory 
approach  is  sometimes  more  physically  intuitive  than 
a  detailed  statistical  analysis,  making  an  interpreta¬ 
tion  of  the  results  more  physically  meaningful.  For 
instance,  it  is  clear  from  this  analysis  that  the  MTF 
describing  reabsorption,  T k(^)>  appears  in  the  cross- 
spectral  density  term  since  light  emitted  remotely  is 
correlated  with  light  emitted  locally  when  a  K  x  ray  is 
produced.  It  appears  in  the  first  power  since  it  appears 
in  only  one  of  two  correlated  paths. 


IV.  CONCLUSIONS 

The  DQE  is  an  important  indicator  of  the  perfor¬ 
mance  of  medical  imaging  systems.  Recent  devel¬ 
opments  in  understanding  noise  transfer  in  medical 
imaging  systems  has  resulted  in  a  generalized  transfer- 
theory  approach  that  can  be  used  to  describe  the  DQE 
and  other  metrics  of  system  performance  for  many 
imaging  systems. 
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TABLE  I.  Values  used  to  determine  r(i/)  from  Metz  and  Vyborny.  Ek  =  59.3  keV. 


Variable 

Value* 

Physical  Meaning 

C 

0.85 

probability  that  an  interacting  x  rayundergoes  a  K-shell  interaction 

u 

0.93 

fluorescent  yield  of  K-shell  photoelectricinteractions 

ip 

0.866 

fraction  of  incident  x-ray  energy  transferredto  K  x  ray,  Ek/Ex 

fK 

0.20 

probability  that  a  K  x  ray  is  reabsorbed  in  thescreen,  depends  on  geometry 

aSee  Ref.  [34] 


Until  recently,  the  transfer- theory  approach  has 
been  limited  to  the  description  of  serial  cascades  of 
quantum  gain  and  quantum  scattering  processes.  As 
part  of  a  program  developing  new  transfer-theory 
relationships,  we  describe  how  parallel  cascades  of 
image-forming  processes  can  be  incorporated  into  the 
transfer-theory  approach.  Parallel  cascades  are  re¬ 
quired  when  more  than  one  image-forming  process 
combines  to  create  the  final  image.  It  has  been  shown 
in  this  article  that  parallel  cascades  can  be  used  with 
the  introduction  of  the  cross  covariance  between  cas¬ 
cades.  A  general  expression  for  the  cross  covariance  of 
correlated  point  processes  has  been  developed,  and  in 
particular,  the  cross  covariance  of  two  amplified  point 
processes  descending  from  randomly  selected  quanta  in 
a  common  input  image  was  examined  which  has  par¬ 
ticular  importance  for  the  analysis  of  medical  imaging 
systems  [Eq.  (81)  for  WSS  conditions].  It  has  been 
shown  how  a  complete  imaging  system  can  be  rep¬ 
resented  in  terms  of  a  schematic  diagram  describing 
relationships  between  elementary  image-forming  pro¬ 
cesses. 

Under  wide-sense  stationary  conditions,  the  Fourier 
transform  of  the  cross  covariance  is  the  cross  spec¬ 
tral  density  function.  With  it,  transfer-theory  models 
can  be  developed  to  describe  the  noise  power  spec¬ 
trum  in  imaging  systems  that  require  the  use  of  par¬ 
allel  cascades  of  image-forming  processes.  The  ex¬ 
ample  of  reabsorption  of  K  x  rays  in  a  radiographic 
screen  was  described  in  this  article.  It  was  shown  that 
the  transfer-theory  approach  gives  the  same  result  ob¬ 
tained  by  Metz  and  Vyborny  using  a  sophisticated  sta¬ 
tistical  analysis  when  the  same  assumptions  are  made. 
Other  examples  may  include:  a)  double-emulsion  film- 
screen  systems  where  light  may  cross  from  one  emul¬ 
sion  to  the  other;  b)  portal  imaging  systems  where 
high-energy  x  rays  may  generate  different  kinds  of  sec¬ 
ondary  quanta  in  the  detector  such  as  electrons  and 
light;  and,  c)  flat-panel  active  matrix  detectors  where 
scattered  light  may  contribute  a  non-negligible  frac¬ 
tion  of  the  image  signal. 

Extension  of  the  transfer-theory  approach  to  include 
parallel  cascades  increases  the  number  of  theoretical 
“tools”  available  to  scientists  and  engineers  in  the 
transfer-theory  “tool-box.”  We  are  currently  devel¬ 
oping  additional  tools  for  the  analysis  of  new  digital 
imaging  systems.  These  new  tools  are  required  to  de¬ 


scribe  noise  transfer  in  the  presence  of  other  situations 
such  as  non-uniform  detector-element  gain  and  sen¬ 
sitivity,  random  errors  in  detector-element  offset  and 
gain  corrections,  polyenergetic  x  rays,  thick  phosphors 
with  variable  conversion  factors  due  to  Swank  noise 
and  the  Lubberts  effect,  and  detectors  with  lag. 
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APPENDIX  A:  ELEMENTARY  PROCESSES 

Medical  x-ray  imaging  system  modeling  is  based  on 
three  elementary  processes.  Under  wide-sense  station¬ 
ary  (WSS)  conditions,  the  elementary  processes,  and 
the  transfer  properties  of  signal  and  noise,  are  summa¬ 
rized  in  this  Appendix. 

1.  Quantum  Gain  and  Selection 

Rabbani  et  al  [13]  described  the  transfer  of  signal 
and  noise  through  a  stochastic  quantum  gain  stage, 
characterized  by  a  random  variable  g  which  has  a  mean 
value  g  and  variance  Barrett  et  al  [19,20]  consid¬ 
ered  this  elementary  process  as  the  limiting  case  of 
amplified  point  process.  They  showed  that  the  uni¬ 
form  mean  value  of  distribution  of  quanta  (i.e.,  mean 
number  of  quanta  per  unit  area)  and  corresponding 
NPS  are  transferred  according  to 

Qout  =  9  Qin  (-^1) 

and 

NPSout  (l v)  =  f  NPSin  (I/)  +  a\qin ,  ( A2) 

respectively.  If  g  is  modeled  as  a  deterministic  gain 
factor,  Eq.  (Al)  remains  the  same,  but  Eq.  (A2)  be¬ 
comes 
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KPSout(v)=g2NPSin(v),  (A3) 

which  lacks  the  second  term  in  Eq.  (A2). 

For  the  special  case  where  g  represents  a  binary  se¬ 
lection  process  such  as  the  responsive  quantum  effi¬ 
ciency  of  the  radiographic  screen,  g  can  have  a  value 
of  0  or  1  only,  0  <  g  <  1,  and  a2  =  g(  1  -  g).  In  this 
case,  it  is  also  possible  to  express  signal  transfer  as 

Qoutir)  =  gqin(r)  (A4) 

where  gin(r)  and  qout(r)  are  sample  functions  of  the 
input  and  output  random  point  processes. 

2.  Quantum  Scatter  (Relocation) 

Quantum  scatter  is  a  translated  point  process  [33] 
whereby  a  quantum  is  randomly  relocated  by  the  dis¬ 
placement  vector  r  and  where  p(r)  is  a  point  spread 
function  (PSF)  describing  the  translation  probability 
density  [13,38]. 

The  transfer  function  through  this  process  have  been 
described  by  both  Rabbani  et  al  [13]  and  Barrett  et 
al  [19,20],  given  by 

Qout  ( r )  =  p(r)  *s  gin  ( r ) ,  ( A5) 

where  *s  represents  the  scattering  process  and  p(r) 
has  unity  area.  The  statistics  of  signal  and  noise  are 
transferred  according  to 

Qout  ~  Qini  (A6) 

and 

NPSoutM  =  [NPS,„(i/)  -  Qin) |X(^)|2  +  qin  (A7) 
where  T(i/)  is  the  Fourier  transform  of  p(r). 

3.  Linear  Filter  (Convolution) 

The  transfer  relationships  through  a  linear  shift- 
invariant  filter  are  described  by  the  convolution  inte¬ 
gral,  given  by  [2] 

d0ut(r)  =  p{r)  *  qin(r),  (A8) 

dout  =  |T(0)|gins  (A9) 

and 

NPS0U(  {v)  =  NPSir»  |T(i/)  |2  (A10) 

where  *  represents  a  convolution,  p(r)  is  the  impulse 
response  commonly  referred  to  as  the  blur  PSF  and 
T(i/)  is  the  characteristic  transfer  function  of  the  fil¬ 
ter,  given  by  the  Fourier  transform  of  p{r).  Unlike 
the  scattering  process,  the  ouput  dout( r)  from  a  linear 
filter  is  not  a  point  process. 


APPENDIX  B:  STATISTICS  OF  PARALLEL 
PROCESSES 

When  two  random  processes  contribute  to  an  out¬ 
put  signal,  the  result  is  a  random  process  that  is  the 
sum  of  two  random  processes.  For  example,  consider 
a  stochastic  system  described  by  c(r)  =  a(r)  4-  b(r). 
The  autocorrelation  of  the  sum  process  is  [39] 

Rc(r,r')=E{c(r)c*(r')} 

=  E{[o(r)+6(r)][o*(r,)  +  6*(r')]} 

=  E{d(r)a*(r')}  +  E{6(r)6*(r')} 

+E{o(r)6*(r')}  +E{&(r)a*(r')} 

=  Ra(r,r')  +Rb{r,r')  +  Rab(r,r')  +  Rba(r,r'),  (Bl) 

where  *  denotes  a  complex  conjugate.  When  a(r)  and 
b(r)  are  both  wide-sense  stationary  (WSS),  c(r)  is  also 
WSS,  and  the  autocorrelation  of  c(r)  in  Eq.  (Bl)  can 
be  written  as 

Rc(t)  =  Ra(r)  +  Rft(r)  +  Raft(r)  +  Rf,a(-r).  (B2) 

The  autocovariance,  K c(r,r'),  is  similar  to  the  auto¬ 
correlation  but  with  the  mean  subtracted.  Similar  to 
above, 

Kc(r,r')  =  E{[a(r)  +  6(r)][o*(r')  +  6*(r')]} 

-E{a(r)  +  5(r)}E{d*(r')  +  6*(r')} 

=  Ka(r,r')  +  K»(r,r')  +  Ka6(r,r')  +  Kfca(r,r'),  (B3) 

and  for  WSS  processes 

Kc(r)  =  K0(r)  +  K»(r)  +  Ko6(r)  +  K6a(r).  (B4) 

The  corresponding  NPS  of  c(r)  is  therefore 
NPS»  =jr{Kc(r)} 

=  NPSa(i/)  +  NPS6(i/)  +  NPSat(i/)  +  NPStaM,  (B5) 

where  each  term  in  Eq.  (B5)  is  the  Fourier  transform 
of  the  corresponding  term  in  Eq.  (B4).  The  terms 
NPS abW  and  NPS6a(i/)  are  cross  spectral  densities 
and  reflect  the  spatial-frequency  dependence  of  the 
auto- covariance. 

By  definition  of  the  cross  covariance  for  two  random 
procesess,  we  have 

Ka6(r,r')=K6*a(r',r).  (B6) 

Under  WSS  conditions,  Eq.  (B6)  becomes 
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Kat(r)  =  K  K-r),  (B7) 

where  r  =  r-r'.  The  cross  spectral  density  NPSafc(i/) 
is  the  Fourier  transform  of  Ka(,(r)  given  by 

NPSo6(r/)  =  ^{Kai,(r)}  =  J  K ab{T)e~^UT dr . 

(B8) 

Substituting  Eq.  (B7)  into  Eq.  (B8),  we  have 
NPSa*(r/)=  [  K la{-r)e~^VT dr 

J  OO 


=  |jF{K6o(r)}|*=NPSLM.  (B9) 

Thus,  we  conclude  that  the  cross  spectral  densities 
NPSat,(i/)  and  NPSf,a(^)  are  conjugate  pairs,  and  the 
sum  of  NPSafc(^)  +  NPSta(^)  is  always  a  real  value, 
i.e., 

NPSo6(i/)  +  NPSfca(i/)  =  2Re{NPSa6(i/)} 

=  2Re{NPS6a(i/)},  (BIO) 

where  Re{  }  denotes  the  real  part  of  a  complex  quan¬ 
tity. 

APPENDIX  C:  STATISTICAL  AVERAGES 
AND  MARGINAL  PROBABILITY  DENSITIES 

Consider  a  continuous  random  variable  y  =  /(x), 
where  x  is  a  continuous  random  variable  with  the  prob¬ 
ability  density  function  Px{x).  The  function  f(x )  is  a 
single- valued  function  of  x,  and  maps  random  variable 
x  into  random  variable  y.  We  define  the  mean  of  the 
continuous  random  variable  y  by  the  equation  [2] 

y  =  E  {y}  =  E|c/(5)}  =  J  g{x)pz{x)dx,  (Cl) 

which  is  the  statistical  average  over  random  variable 
x.  If  the  probability  density  function  pg(y)  is  defined, 
then 

y  =  E{y}  =  [  yPy(y)dy.  (C2) 

Eqs.  (Cl)  and  (C2)  give  the  same  result  for  the  statis¬ 
tical  average  of  y  —  g(x).  Equation  (Cl)  refers  to  the 
sample  space  of  x  and  Eq.  (C2)  refers  to  the  sample 
space  of  y . 


For  a  set  of  several  random  variables,  the  statistics 
of  each  are  called  marginal  [2].  For  example,  Px(x) 
is  the  marginal  probability  density  of  x  and  py  (y)  the 
marginal  probability  density  of  y  for  their  joint  den¬ 
sity  function  Px,y{x,  y).  The  following  relationships  are 


satisfied  [2]  for  pz(x)  and  Py(y), 

Pi(x)  =  /  Pi,y(x,y)dy 

J  — OO 

(C3) 

and 

Py(y)=  /  Px,y(x,y)dx> 

J  — OO 

(C4) 

which  are  the  property  of  marginal  densities. 
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Appendix  II:  Scatter  Operator 

This  appendix  consists  of  a  manuscript  completed  this  year  describing  how  x-ray  scatter  can  be 
incorporated  into  the  linear-systems  approach  as  a  “scatter  operator.”  It  has  been  submitted  to  Medical 
Physics  for  publication. 
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ABSTRACT 

Rabbani,  Shaw  and  Van  Metter  described  the  transfer  of  the  auto-covariance  and  image  noise-power 
spectrum  through  a  photon  scattering  (image  blurring)  process.  We  incorporate  their  result  into  a  linear- 
systems  transfer-theory  framework  by  describing  a  generalized  scatter  operator  acting  with  an  associated 
point-spread  function.  Scatter  is  a  translated  point-process,  and  it  is  shown  that  this  operator  is  not 
distributive  over  addition  or  multiplication.  While  in  general  it  does  not  commute  with  addition  or  mul¬ 
tiplication,  it  does  commute  with  itself  and  with  a  binomial  selection  process  (as  might  represent  the 
quantum  efficiency  of  a  detector).  Strictly  speaking  it  is  not  associative,  although  a  similar  property  that 
is  true  is  described.  This  formalism  is  used  to  explain  the  degradation  of  the  detective  quantum  efficiency 
(DQE)  caused  by  both  the  x-ray  focal- spot  size  and  optical  blur  in  a  radiographic  screen  as  illustrative 
examples.  It  is  shown  that  for  quantum- noise  limited  systems  in  general,  the  frequency  dependence  of 
the  system  DQE  is  a  direct  consequence  of  the  statistical  properties  of  image  quanta  as  represented  by 
the  scatter  operator.  This  extended  linear-systems  approach  has  wide  application  for  the  description  of 
imaging  systems  in  general. 

Keywords:  linear-systems,  transfer  theory,  detective  quantum  efficiency,  DQE,  convolution,  modulation- 
transfer  function,  MTF,  noise-power  spectrum,  NPS,  cascaded  systems 

I.  INTRODUCTION 

The  design  and  analysis  of  medical  imaging  systems  requires  a  mathematical  understanding  of  important 
principles  and  processes  that  contribute  to  image  formation.  Linear-systems  transfer  theory  is  an  approach 
used  by  many  authors  to  describe  the  performance  of  linear  and  shift-invariant  (LSI)  imaging  systems  in 
both  the  spatial  and  spatial- frequency  domains.1"4  Image-blurring  mechanisms  are  described  in  terms 
of  linear  filters1’5,6  as  a  convolution  with  a  specified  point-spread  function  (PSF)  in  the  spatial  domain, 
or  alternatively  by  multiplication  with  the  modulation-transfer  function  (MTF)  in  the  spatial- frequency 
domain.7,8  Complex  systems  can  be  represented  as  a  serial  cascade  of  simple  processes,  or  stages,  where 
the  output  from  one  stage  forms  a  virtual  input  to  the  next.  Excellent  early  review  papers  on  this  topic 
have  been  written  by  Doi,  Rossmann  and  Haus,9  and  Metz  and  Doi.2 

This  linear-filter  approach  describes  only  the  expectation  (noise-free)  effect  of  image- blurring  mecha¬ 
nisms  since  it  represents  blur  as  a  deterministic  process.  In  reality,  blur  generally  results  from  the  random 
scattering  of  individual  image  quanta,  and  the  effect  of  this  scatter  on  image  noise  must  be  described 
separately  using  stochastic-process  theories.2’10  For  the  special  case  of  wide-sense-stationary  (WSS)  noise 
processes,11  this  noise  can  be  characterized  in  terms  of  the  auto-covariance,  or  its  Fourier  transform,  the 
noise-power  spectrum  (NPS).  Dainty  and  Shaw,8  Barrett  and  Swindell,1  Macovski,  Parker  and  others 
have  used  this  approach  extensively  in  their  texts  to  describe  important  principles  and  fundamental  prop¬ 
erties  of  many  system  designs. 
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The  distinction  between  scatter  and  deterministic  blur  has  been  recognized  by  many  investigators, 
including  Dainty  and  Shaw,8  Wagner,12  Metz  and  Doi,2  Sandrik  and  Wagner,13  Metz  and  Vyborny,14 
and  Barrett  and  Swindell.1  The  NPS  of  a  random  process  undergoing  a  deterministic  blur  is  “passed 
through”  the  square  of  the  blur  MTF  (e.g.  Barrett  and  Swindell1  Eq.  3.200,  Parker6  Eq.  24.28).  More 
recently,  transfer  of  the  NPS  through  a  scattering  process  was  described  by  Shaw  and  Van  Metter,15 
derived  theoretically  by  Rabbani  et  al.16  using  moment-generating  functions,  and  later  by  Barrett  et  al.17 
using  point-process  theory.  They  note  that  the  correlated  component  of  the  NPS  is  passed  through  the 
squared  scatter  MTF  while  the  uncorrelated  component  is  not.15-16  Their  result  has  been  experimentally 
demonstrated  by  Maidment  and  Yaffe,18  and  derived  in  a  different  form  by  Mulder.19 

The  scatter  result16  has  been  used  by  Nishikawa  et  al.20’21  and  others  in  a  small  number  of  applications, 
but  its  implications  are  of  a  more  general  nature.  We  incorporate  scatter  into  the  linear-systems  approach 
by  describing  a  generalized  scatter  operator  as  an  alternative  to  the  convolution  operation  that  is  used 
to  describe  deterministic  blur.  The  word  “convolution”  means  to  fold  or  roll  together  (Webster),  and  is 
used  to  represent  the  blur  of  a  linear  filter  when  the  filter  PSF  is  “folded  into”  the  image  signal.  In  some 
sense,  scatter  can  be  interpreted  as  a  “stochastic  convolution,”  as  this  folding  together  is  a  consequence 
of  the  random  relocating  of  individual  image  quanta.  Use  of  the  scatter  operator  (and  hence  the  Rabbani 
NPS)  as  a  substitute  for  the  convolution  integral  to  describe  image  blurring  processes  means  that  linear- 
systems  transfer  theory  can  be  used  to  describe  both  image-signal  and  image-noise  (second  order)  transfer 
characteristics  in  a  cascaded-systems  analysis. 

Under  certain  circumstances,  the  difference  between  scatter  and  deterministic  blurring  is  insignificant. 
For  instance,  Barrett  and  Swindell22  note  that  when  x  rays  are  converted  to  light  in  a  radiographic  screen, 
the  statistical  nature  of  the  light  photons  can  be  ignored  if  a  large  number  of  photons  are  generated 
and  hence  each  x  ray  produces  light  that  is  adequately  described  by  the  deterministic  PSF.  The  work 
described  here  discusses  this  in  more  detail,  and  extends  this  observation  to  show  that  it  is  spatial-frequency 
dependent.  In  addition,  if  the  light  is  subsequently  collected  by  an  inefficient  optical  system,  the  statistical 
nature  of  the  optical  photons  cannot  be  ignored.  This  will  occur,  for  example,  if  only  a  few  light  quanta 
are  detected  per  interacting  x  ray  since  it  is  not  possible  to  know  where,  within  each  PSF  distribution, 
the  detected  quanta  originated.  Thus,  when  assessing  new  systems  or  system  designs,  it  is  necessary 
to  correctly  account  for  these  statistical  properties.  Ignoring  them  can  result  in  much  wasted  effort  in 
designing  new  systems  that  have  no  chance  of  achieving  a  sufficiently  high  detective  quantum  efficiency 
(DQE). 


II.  BACKGROUND 

Linear  systems  can  be  described  in  terms  of  the  transfer  of  image  signal  and  noise  from  input  to  output. 
In  this  context,  an  image  is  a  spatial  distribution  of  quanta  (e.g.  x  rays  or  light).  Signal  transfer  corresponds 
to  the  transfer  of  the  expected  distribution  of  image  quanta.  Noise  transfer  corresponds  to  the  transfer  of 
parameters  that  describe  image  noise,  such  as  the  autocovariance  or  NPS.  In  this  article,  a  third  type  of 
transfer  relationship  will  be  used  that  we  call  the  sample  transfer,  describing  the  transfer  of  a  particular 
sample  distribution  of  image  quanta. 

Quanta  have  negligible  spatial  extent,  and  one  way  of  describing  a  sample  distribution  of  quanta  is 
as  q{r),  a  superposition  of  many  ^-functions  where  each  5-function  represents  one  quantum.  While  this 
approach  is  often  convenient,  5-functions  are  generalized  functions3,23  and  must  be  treated  with  care.  In 
addition,  it  is  convenient  to  write  the  ensemble  average,  or  expectation,  of  q(r)  as  E{g(r)}.  Thus,  while 
q(r)  is  a  sample  distribution  of  quanta  representing  a  particular  image,  ¥j{q(r )}  describes  the  expectation 
number  of  quanta  per  unit  area  as  a  function  of  position,  and  can  be  viewed  as  a  noise-free  image.  Both 
have  dimension  of  quanta  per  unit  area. 

The  generalized  function  q{r)  is  not  directly  measurable.  Measurable  quantities,  such  as  the  output 
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quanta  which  light  quanta  quanta  quanta  which  are  elements 
interact  in  the  detected 

screen 


Figure  1:  The  hypothetical  system  represented  as  a  cascade  of  linear  stages. 


from  a  detector  element  (see  Appendix  I)  are  expressed  as  d  given  by 


where  Ad  is  the  area  of  the  detector  element  and  g0  is  a  conversion  factor  relating  the  number  of  interacting 
quanta  to  the  detector  signal.  When  Ad  is  expressed  as  Ad(r),  the  area  of  a  hypothetical  element  centered 
at  position  r,  the  signal  d  is  written  as  d(r)  and  has  been  called  the  “presampling”  detector  signal  as  it 
is  a  function  that  describes  the  detector  signals  when  evaluated  (“sampled”)  at  positions  corresponding 
to  the  centers  of  the  detector  elements.  In  the  following,  we  use  q  to  represent  a  distribution  of  quanta 
and  d  to  represent  a  physically  measured  -  or  measurable  -  detector  signal  in  terms  of  an  analog-to-digital 
converter  (ADC)  value.  Image  noise  is  defined  in  terms  of  the  noise-power  spectrum  (NPS)  in  q(r)  for  the 
special  case  of  E{q(r)}  equal  to  a  constant  throughout  the  image  to  satisfy  the  requirement  for  wide-sense 
stationary  noise  processes  where  the  mean  and  autocorrelation  are  invariant  with  r.11 

The  utility  of  the  linear-systems  approach  is  illustrated  by  analyzing  the  simplified  hypothetical  system 
shown  in  Fig.  (1)  to  determine  the  presampling  detector  signal.  The  system  consists  of  a  thin  radiographic 
screen  coupled  to  a  CCD  camera  through  a  lens  assembly.  All  components  are  assumed  to  be  linear 
and  shift  invariant.  The  screen  quantum  efficiency  is  a  (assumed  to  be  unity  for  the  moment)  and  each 
interacting  x  ray  generates  exactly  m  optical  quanta  at  some  point  of  interaction  in  the  screen.  The  optical 
quanta  are  scattered  (i.e.  spread)  before  leaving  the  screen  according  to  a  PSF  Pj.(t’).  The  optical  quanta 
are  coupled  to  the  CCD  camera  with  a  lens  having  a  collection  efficiency  B.  The  camera  integrates  incident 
optical  quanta  in  discrete  detector  elements  and  the  signals  from  all  elements  are  combined  to  generate  a 
digital  image.  This  model  does  not  describe  the  effect  of  such  things  as  variable  x-ray  energy,  or  variable 
interaction  depth  or  gain  in  the  screen.  As  such,  it  is  a  simplified  model  and  is  chosen  primarily  to  highlight 
the  difference  between  use  of  the  conventional  and  “stochastic”  convolutions  in  linear-systems  theory. 

The  serial  nature  of  the  model  is  illustrated  in  Fig.  (1)  where  q(r)  represents  the  distribution  of  x-ray 
quanta  incident  on  the  screen.  The  distribution  of  interacting  quanta  is  qi{r)  and  the  distribution  of  optical 
quanta  as  generated  in  the  screen  is  <?m(r)  where 

E{tfm(r)}  =  E  {qi{r)}m  =  E  {q{r)}am.  (2) 

The  two-dimensional  NPS  in  the  distribution  of  optical  quanta  is  NPSm(fc)  given  by16 


NPSm(fc)  =  m2NPS/(fc)  =  m2qi 


(3) 
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where  k  is  a  two-dimensional  frequency  vector  (with  units  cycles  per  unit  length)  in  the  direction  of  r 
and  the  NPS  of  qj{r)  is3  NPSj(fc)  =  ?/•  The  assumption  of  a  deterministic  gain  m  is  not  always  a  valid 
assumption,  and  is  discussed  below. 

If  it  can  further  be  assumed  that  each  interacting  x-ray  quantum  produces  exactly  the  same  distribution 
of  light  at  the  screen  surface,  and  that  no  optical  quanta  are  lost  or  absorbed  in  the  screen,  then  the 
expectation  distribution  of  g0(r),  the  optical  quanta  that  exit  from  the  screen,  can  be  expressed  using 
linear-systems  theory  as 

/°° 

E{^m(r')}Pr(^  -r')dV  (4) 

-00 

where  **  represents  a  two-dimensional  convolution  integral.3  A  pure  scattering  process  neither  creates 
nor  destroys  quanta,  and  hence  [),.  (r)  must  necessarily  be  normalized  to  unity  area.  Thus,  as  shown  in 
Eq.  (4),  the  expectation  blurred  optical  image  E{y0(r)}  can  be  expressed  as  a  convolution  of  an  expectation 
“pre-blurring”  optical  image  E{<?m(r)}  with  the  appropriate  PSF. 

The  relationship  described  by  Eq.  (4)  can  be  expressed  in  the  spatial-frequency  domain  as 

E{Q0(k)}  =  E{QI(k)}mTr(k)  (5) 

where  E{<30(fc)},  E{Q[{k)}  and  T r(k)  are  the  two-dimensional  Fourier  transforms  of  E{y0(r)},  E{g/(r)} 
and  pr  (r)  respectively.  The  two-dimensional  MTF  for  this  system  is  therefore  given  by 

MTF(fc)  =  |Tr(fc)|  (6) 


which  has  a  value  of  unity  when  k  =  0  since  pr  ( r )  is  required  to  have  unity  area. 

The  next  stage  in  the  model  is  propagation  of  the  optical  quanta  and  integration  in  the  digital  detector. 
In  Appendix  I  it  is  shown  that  the  presampling  detector  signal  d(r)  can  be  represented  as  a  convolution 
of  the  distribution  of  interacting  quanta  with  an  appropriate  rectangular  aperture  function  having  unity 
height.  An  overall  expression  for  the  expectation  presampling  detector  signal  in  Cartesian  coordinates  can 
therefore  be  written  as 


E{d(r)}  =  E{d(:r,y)}  =  g0 


[E  {qi(x,  y)}m  *  *  pr(z,  y)]  /3  *  *  II 


(7) 


where  f3  describes  the  fraction  of  optical  quanta  leaving  the  screen  that  are  integrated  in  the  detector,  ax 
and  ay  are  the  x-  and  y-direction  dimensions  of  the  active  regions  of  individual  detector  elements. 

The  convolution  operator  commutes  with  both  itself  and  multiplication,3  and  hence 


E {d{x,y)}  =  g0mPE{qi(x,y)}  *  * 


pr(x,y)**II 


*jJL\ 

(Lx  &y  J 


(8) 


This  relationship  can  be  written  in  the  spatial-frequency  domain  as 

E{D{u,v)}  =  g0m/3axay  E{Qi(u,  v)}  Tr(it,u)  sinc(7raxu)  sinc(7r ayv)  (9) 

where  E{  D(u,  v)}  is  the  (complex)  expected  value  of  the  Fourier  transform  of  d(x,  y),  u  and  v  are  spatial  fre¬ 
quencies  in  the  x-  and  y-directions  respectively,  and  sinc(0)  =  sin(0)/ 6 .  The  product  y o o-x °j/ s i n c ( ti' ) s i n c (  ,t a y 
has  been  called  the  sampling  aperture  optical  transfer  function  by  Giger  and  Doi.24  The  two-dimensional 
MTF  for  this  system  including  the  effect  of  integration  in  detector  elements  is  therefore  given  by 

MTF(u,u)  =  |Tr(u,u)  sinc(7raxu)  sinc^a^)!.  (10) 


While  x-ray  quanta  are  always  Poisson  distributed  when  distributed  with  a  uniform  mean,  this  is  not 
generally  so  with  other  image  quanta  such  as  light  from  a  screen.  Therefore,  Eq.  (7)  cannot  be  used  directly 
to  determine  the  NPS  in  d(x,y).  This  is  in  part  because  the  convolution  integral  correctly  describes  the 
transfer  of  mean  values,  but  not  noise,  through  a  stochastic  system.  In  the  following  sections,  the  “stochastic 
convolution”  representation  of  scatter  is  suggested  as  a  means  of  bringing  the  Rabbani  result  to  this  type 
of  analysis. 
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III.  THEORY 

The  blurring  of  an  image  when  described  by  the  convolution  integral  is  given  by  Eq.  (4)  for  LSI  systems. 
In  this  expression,  pr(r  -  r')  describes  the  fraction  of  the  signal  in  an  infinitesimal  element  d2r  that  is 
“mis-located”  to  the  new  position  r  as  a  result  of  the  optical  blurring  process.  The  NPS  is  transferred 
through  a  convolution  integral  according  to3 

NPSd(fc)  =  NPS(fc)  |T(fe)|2  (11) 

where  NPS(fc)  is  the  initial  NPS  and  T (k)  is  the  characteristic  transfer  function  of  the  blurring  process. 
We  will  call  this  process  “deterministic  blur,”  as  indicted  by  the  subscript  d,  to  emphasize  the  deterministic 
nature  of  the  convolution.  However,  this  result  is  not  generally  used  to  describe  the  NPS  transfer  through 
imaging  systems  because  it  does  not  accommodate  the  discrete  nature  and  statistical  properties  of  the 
image  quanta  responsible  for  blurring  processes. 

Rabbani  et  al.16  developed  a  theoretical  description  of  the  resulting  NPS  when  independent  quanta  are 
scattered  according  to  a  probability  density  distribution  pr(r)  that  specifically  accounts  for  the  statistical 
properties  of  these  quanta.  When  pr(r)  is  shift  invariant,  a  general  expression  for  the  resulting  NPS 
NPSs(fc)  is  given  by 

NPSj(fc)  =  [NPS (fc)  -  q]  |T(fc)|2  +  q  (12) 

where  q  is  the  average  distribution  of  image  quanta  before  the  blur.  We  will  call  this  “stochastic  blur^  as 
indicated  by  the  subscript  s.  This  result  shows  that  the  correlated  component  of  the  NPS,  NPS(fe)  -  q,  is 
passed  though  the  squared  MTF  while  the  uncorrelated  component,  q,  is  not.16  A  correlated  component 
can  result  from  a  conversion  with  gain  from  one  form  of  quanta  to  another  as  described  in  Sec.  III.B.  The 
NPS  resulting  from  stochastic  blur  can  also  be  expressed  in  terms  of  the  NPS  resulting  from  deterministic 
blur  [Eq.  (11)]  as 

NPSs(fc)  =  NPS d(k)  +  [l  -  |T(fc)|2]  q  (13) 

showing  that  at  low  frequencies  where  |T(fc)|2  »  1,  NPSd(fc)  and  NPS,(fc)  are  similar  and  hence  the 
deterministic  and  stochastic  results  are  similar.  They  are  also  similar  if  the  distribution  of  input  quanta 
are  statistically  correlated,  implying  that  NPS(fc)  >  q  (Sec.  III.B),  except  at  very  high  spatial  frequencies. 
However,  at  non-zero  frequencies  where  |T(fc)l  <  1,  the  deterministic  result  in  general  underestimates 

image  noise. 

Use  of  linear-systems  theory  in  general,  and  the  convolution  operator  in  particular,  to  describe  image 
blurring  processes  is  a  very  powerful  technique.  Thus,  while  the  convolution  integral  results  in  an  NPS 
described  by  Eq.  (11),  it  seems  appropriate  to  introduce  a  modified  operator  that  will  result  in  the  NPS 
described  by  Eq.  (12).  It  is  intended  that  this  modified  operator  unify  the  description  of  both  image 
blur  and  image  noise  within  a  linear-systems  framework.  That  is,  relatively  simple  expressions  similar  to 
Eq.  (7)  can  be  developed  to  describe  both  image-signal  and  image- noise  transfer  characteristics  of  particular 
systems  by  using  this  scatter  operator. 

A.  The  Scatter  Operator,  *s 

The  input  to  an  imaging  system  is  a  distribution  of  quanta,  q(r ),  which  may  be  written  as 

«(r)  =  E  *(’—•<)  ^ 

*= l 

where  N  is  a  Poisson  random  variable  describing  the  total  number  of  incident  quanta  and  r*  is  a  vector 
describing  the  position  of  the  ith  quantum.  Thus,  q(r)  is  a  sample  distribution  of  a  random  process 
generating  image  quanta  with  the  appropriate  statistical  distribution. 
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Convolution  Scatter 


Figure  2:  Deterministic  blur  represented  with  a  convolution  integral,  and  scatter,  are  shown  in  the  left 
and  right  columns  respectively.  The  convolution  “smears”  the  input  according  to  the  PSF  (a  to  c),  while 
scatter  “relocates”  the  input  quanta  according  to  the  probability  distribution  given  by  the  PSF  (d  to  f). 

Rabbani  et  al.  derived  Eq.  (12)  in  terms  of  the  set  of  transition  probabilities  {p^},  where  is  the 
probability  that  a  quantum  in  pixel  j  is  scattered  to  pixel  k.  For  LSI  systems  this  transition  probability  is 
independent  of  position  and  hence  the  probability  of  being  scattered  to  a  pixel  displaced  by  the  vector  Ar 
can  be  written  pAr-  The  distribution  of  quanta  following  a  scattering  process  is  described  by  qs(r)  where 

N 

qs(r)  =  ^2  S  (r  -  [r j  +  Ar*])  (15) 

i— 1 

and  where  Ar*  describes  the  scatter  displacement  of  the  ith  quantum.  The  displacement  A r i  is  a  random 
vector  variable  with  a  set  of  probabilities  {pat**}  that  describe  the  scatter  PSF  p(r)  where 

PAri  =p(r)d2r|r=Ari  (16) 

and  where  d2r  is  the  area  of  an  infinitesimal  image  element.  Equation  (15)  describes  the  “stochastic 
convolution”  of  q(v)  with  p(r),  and  is  written  in  short  form  for  two-dimensional  scatter  as 

qs{r)  =q{r)*s*sP(r).  (17) 

The  scatter  operator  describes  transfer  of  a  sample  distribution  of  image  quanta.  It  is  a  translated  point- 
process,  similar  to  the  translated  Poisson-process  described  by  Snyder  and  Miller,10  although  not  restricted 
to  Poisson  processes. 

Thus,  while  the  convolution  integral  “smears”  the  input  according  to  the  PSF,  scatter  “mis-locates 
input  quanta  according  to  a  probability  described  by  the  PSF.  This  difference  is  illustrated  in  Fig.  (2). 
In  the  column  on  the  left,  the  input  (2a)  consists  of  randomly  positioned  quanta,  each  represented  as  a 
<5-function.  This  input  is  convolved  with  the  PSF  (2b)  producing  the  result  (2c)  with  a  corresponding  NPS 
described  by  Eq.  (11).  On  the  right,  the  same  input  undergoes  a  scatter  resulting  in  a  random  relocation 
of  all  quanta  (2f)  according  to  the  same  PSF.  The  NPS  of  this  result  is  described  by  Eq.  (12). 

1.  Conditions  of  Use 

The  conditions  of  use  of  the  scatter  operator  *5  are  essentially  the  assumptions  described  by  Rabbani 
et  al.  in  the  derivation  of  Eq.  (12).  They  include: 
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1.  q(r)  must  represent  a  distribution  of  independent  quanta,  with  units  of  quanta  per  unit  area  (a 
one-dimensional  form  would  have  units  of  quanta  per  unit  length); 

2.  the  scattering  process  is  shift  invariant; 

3.  all  scattering  events  are  independent; 

4.  p (r)  represents  the  probability  (per  unit  area)  that  any  given  quantum  is  scattered  (relocated)  by  a 
displacement  vector  r;  and, 

5.  quanta  can  be  neither  created  nor  destroyed  by  the  scattering  process  [the  integral  of  p(r)  over  all 
space  is  unity]. 

2.  Numerical  Implementation 

It  is  intended  that  *s  be  used  symbolically  [e.g.  Eq.  (17)]  to  represent  scattering  processes,  and  the 
Rabbani  equation  be  used  to  describe  NPS  transfer.  However,  a  numerical  implementation  is  also  possible 
using  simple  Monte  Carlo  methods25  should  a  numerical  simulation  of  the  blurring  process  be  desired.  The 
steps  required  in  such  a  calculation  are  summarized  below  and  must  be  repeated  for  each  quantum  being 
scattered: 

1.  Choose  the  scattering  polar  angle  6  where  0  <  6  <  27T.  This  will  be  a  uniformly  distributed  random 
variable  if  the  scattering  process  is  isotropic.  This  step  is  not  required  for  a  one-dimensional  geometry. 

2.  Choose  the  scatter  displacement  magnitude  |Arj|  which  is  a  random  variable  with  a  probability 
distribution  given  by  PArt  [Eq.  (16)]. 

3.  Reassign  the  location  of  the  quantum  being  scattered  by  the  displacement  Ar,  determined  from  the 
direction  and  distance  calculated  above. 

Computation  time  is  therefore  dependent  on  the  total  number  of  quanta  in  the  image  at  the  scattering 
stage,  but  not  on  the  size  or  number  of  pixels.  Each  image  shown  in  the  results  section  of  this  paper  contain 
approximately  107  quanta  in  a  256  x  512-pixel  matrix,  requiring  approximately  1  minute  to  implement 
the  scatter  on  a  Sun  SPARC  10. 

B.  Scatter  in  a  Cascaded  Model 

The  operator  *s  represents  a  single  scattering  process.  However,  most  real  systems  are  better  repre¬ 
sented  as  a  cascade  of  multiple  scattering  processes  intermixed  with  quantum  gain  (conversion)  processes 
as  illustrated  in  Fig.  (1).  A  consequence  of  the  gain  stages  is  that  they  can  introduce  spatial  correlations 
into  the  distribution  of  image  quanta  and  thereby  influence  how  image  noise  is  passed  through  scattering 
processes.  This  was  described  by  Rabbani  et  al.,16  who  showed  that  an  increase  in  the  number  of  quanta 
representing  the  image  by  a  gain  stage  having  a  mean  g  and  variance  a2  affects  the  NPS  by 

NPSs(fc)  =  <?2NPS(fc)  +  o2gq.  (18) 

where  NPS(fc)  and  q  are  the  NPS  and  mean  respectively  of  the  input  distribution  of  image  quanta  q(r). 
When  the  distribution  of  quanta  after  the  gain  is  uncorrelated,  the  NPS  is  equal  to  the  mean,  i.e.  NPS9  (k)  = 
gNPS(k)  =  gq.  It  is  clear  that  the  actual  NPS  given  by  Eq.  (18)  is  always  greater  than  this  value  for  g  >  1 
and  non-deterministic  gains  (<r2  >  0).  Thus,  the  NPS  of  the  distribution  of  quanta  following  a  gain  stage 
(in  which  the  quanta  are  “clumped”  according  to  the  distribution  of  quanta  before  the  gain)  is  always 
greater  than  the  NPS  of  the  same  number  of  uncorrelated  quanta.  Equation  (18)  is  an  extension  into  the 
spatial-frequency  domain  of  the  effect  of  a  stochastic  gain  on  the  signal  variance.1,26,27  The  deterministic 
gain  used  in  Eq.  (3)  corresponds  to  a2  =  0. 
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Figure  3:  The  statistical  nature  of  light  photons  generated  in  a  screen  become  less  important,  and  stochastic 
and  deterministic  blurs  merge,  as  the  number  of  photons  generated  per  x  ray  increases:  (a)  1,  (b)  2,  (c) 
10,  (d)  100  and  (e)  1000. 

1.  Cascading  Gain  and  Scatter  -  a  Radiographic  Screen 

In  the  radiographic  screen  example,  each  interacting  x  ray  produces  many  light  photons  that  are 
scattered  before  leaving  the  screen  on  the  exit  side.  If  the  gain  is  assumed  to  be  deterministic  for  simplicity 
with  a  value  5,  the  resulting  NPS,  NPSp?s(fc),  is  obtained  by  combining  Eqs.  (18)  and  (12): 

NPSfl,s(fc)  =  NPS^(fc)  +  gq[  1  -  |T(fc)|2]  (19) 

where  NPSfl  <*(fc)  is  given  by 

NPS9)d(fc)  =  <?2NPS(fc)|T(fc)|2  (20) 

and  is  the  NPS  that  would  be  obtained  if  the  blur  were  deterministic,  obtained  by  combining  Eqs.  (18) 
and  (11). 

Equation  (19)  shows  that  in  general  NPS9>s(fc)  is  greater  than  NPS g,d{k)-  They  are  similar  only 
when  the  first  term  on  the  right  side  of  Eq.  (19)  dominates  which  will  occur  when  pNPS(fc)|T(fe)|2  » 
<7[1  —  |T(fc)|2].  This  may  occur  when  g  is  sufficiently  large  or  when  |T(fc)|2  ~  1  (i.e.  \k\  ~  0).  At  sufficiently 
high  spatial  frequencies  for  which  |T(fc)|2  q/[gHPS(k)  +  <7]  (a  condition  often  true  since  g  can  be  very 
large  and  NPS(fc)  can  be  no  less  than  q),  the  second  term  will  dominate  and  the  Rabbani  NPS  resulting 
from  scatter  is  much  greater  than  the  NPS  resulting  from  deterministic  blur. 

The  above  condition  for  large  gain  factors  is  illustrated  in  Fig.  (3)  where  the  distribution  of  light 
photons  is  calculated  following  a  scatter  and  various  gain  values.  Figure  (3a)  shows  the  same  distribution 
of  interacting  quanta  used  in  Fig.  (2).  The  distribution  of  light  photons  following  scatter  (implemented 
using  the  Monte  Carlo  method)  approaches  the  deterministic  convolution  result  [Fig.  (2c)]  as  g  is  increased 
from  2  [Fig.  (3b)]  to  1000  [Fig.  (3e)].  This  illustrates  how  similar  the  stochastic  and  deterministic  blur 
results  are  when  following  a  large  gain  stage  as  predicted  by  Eqs.  (19)  and  (20). 

2.  Error  Associated  with  the  Deterministic-Blur  Model 

Equation  (19)  shows  that  the  NPS  of  a  system  consisting  of  a  noise-free  gain  followed  by  a  blurring 
process  is  underestimated  when  the  blur  is  represented  as  a  convolution  integral  rather  than  as  a  scatter. 
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MTF  Value 


Figure  4:  The  relative  error  in  the  NPS  following  an  amplification  stage  (such  as  generation  of  light  photons 
in  a  radiographic  screen)  resulting  from  use  of  the  conventional  convolution  for  a  Poisson-distributed  input 
is  shown  as  a  function  of  the  MTF  value  for  gains  of  1,  10,  100  and  200. 


This  relative  error  is  given  by  the  fraction  E(k)  where 


E(k)  = 


NPSfl,s(fc)-NPSg,d(fc) 
NPS  g,d(k) 


L|T(fc)|2 


(21) 


for  an  uncorrelated  Poisson-distributed  input.  Figure  (4)  shows  this  error  increases  dramatically  as  the 
MTF  value  decreases.  In  the  absence  of  gain  ( g  =  1),  the  error  is  approximately  56%  when  the  MTF 
value  is  0.8,  and  increases  to  a  factor  of  99  when  the  MTF  value  is  0.1.  In  a  radiographic  screen,  a  typical 
gain  of  g  =  200  (the  number  of  light  quanta  that  exit  the  screen  per  interacting  x  ray)  reduces  this  error 
to  0.5%  and  46%  respectively.  In  general,  the  error  caused  by  use  of  the  convolution  integral  depends  on 
gains  throughout  the  system,  and  hence  cannot  be  determined  without  considering  the  entire  system. 


C.  Properties  of  the  Scatter  Operator 
1.  Mean  Value 

The  expectation  result  of  scatter  approaches  the  convolution  integral  (noise-free)  result  described  by 
Eq.  (4).  This  was  shown  by  Rabbani  et  al.  in  Ref  16,  Eq.  (22): 

Vj  ==  Y  XkPj-k  (22) 

k 

indicating  that  the  average  number  of  quanta  in  pixel  j  after  the  scattering  process,  yj,  is  equal  to  the 
average  number  of  quanta  in  pixel  k  before  scattering,  multiplied  by  the  transition  probability  Pj-k 
and  summed  over  all  pixels.  Combining  this  with  Eq.  (16)  and  taking  the  limit  of  infinitesimally  small 
pixels  results  in 

/oo 

E{<?(r')}p(r-r')d2r  (23) 

-OO 

which  is  recognized  as  being  the  two-dimensional  convolution  integral  of  the  incident  expectation  dis¬ 
tribution  with  the  scatter  PSF.  This  result  indicates  that  the  linear-filter  approach  describes  the  true 
mean  (noise-free)  image.  The  convolution  integral  may  be  used  to  describe  the  expectation  distribution  of 
quanta  after  a  scatter  process,  but  details  of  the  actual  sample  distribution  are  dependent  on  the  statistical 
properties  of  the  quanta  as  represented  by  the  scatter  operator. 
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2.  Other  Properties 

The  following  results  are  expressed  using  one-dimensional  geometry.  Corresponding  results  are  obtained 
using  two-dimensional  geometry.  In  each  case,  q(x)  represents  an  actual  sample  distribution  of  quanta. 

1.  Unlike  the  convolution  integral,  the  two  terms  do  not  commute,  i.e. 

q(x)  *s  p(z)  ±  p(x)  *s  q{x)  (24) 

although  this  is  not  a  particularly  meaningful  observation  as  one  term  must  be  a  distribution  of 
quanta  and  the  other  a  probability  density  distribution  (a  PSF),  and  the  roles  cannot  be  reversed. 

2.  In  general,  scatter  does  not  commute  with  multiplication,  i.e. 

[gq(x)}  *s  p(z)  ^  9  [q(x)  *s  p(^)]  (25) 

since  the  NPS  of  each  side  evaluated  using  Eq.  (18)  differ.  An  exception  to  this  occurs  when  g 
represents  a  binomial  selection  process  such  as  a  quantum  efficiency.  In  this  special  case,  g  is  a 
random  variable  having  a  value  1  or  0  only  where  the  value  1  occurs  with  a  probability  g,  0  <  g  <  1, 
and  <r2  =  <?(l-g).  The  NPS  of  each  side  are  then  equal,  and  scatter  does  commute  with  multiplication. 
This  means  that  the  order  of  sequential  binomial  selection  and  scattering  processes  can  be  reversed 
without  affecting  the  NPS  of  the  outcome. 

3.  The  scatter  operator  is  not  distributive  under  addition,  i.e. 

[q(x)  *s  Px(x)]  +  [q{x)  *s  p2(x)]  /  Q(x)  *s  [Pi(®)  +  P2(a:)]  (26) 

since  the  PSF  must  have  unity  area  and  it  is  not  possible  for  each  of  Pi  (x)  and  p2(®)  as  well  as 
the  sum  p^(x)  4-  p2 (:r )  to  be  so  normalized.  Nor  is  it  distributive  under  multiplication  for  the  same 
reason. 

4.  Scatter  is  not  associative,  i.e. 

[<j(:r)  *«  Pi  (a)]  *s  p2(z)  ^  Q{x)  *s  [Pifa)  P2(a:)]  (27) 

since  it  is  not  meaningful  to  write  Pi(x)  P2(a:)  35  one  term  must  be  a  distribution  of  quanta. 

A  consequence  of  the  above  results  is  that  it  is  generally  not  possible  to  re-group  or  re-arrange  ex¬ 
pressions  that  include  a  scatter  operation  in  the  same  way  that  it  is  possible  to  re-arrange  expressions 
involving  a  convolution.  This  is  because  transfer  expressions  for  the  NPS  can  be  cascaded,  but  are 
generally  not  multiplicative. 

5.  An  important  property  that  does  apply  to  the  scatter  operator  is 

[<?(z)  *s  Pi  (sc)]  *s  p2(x)  =  q(x)  *s  [Pi(z)  *  P2(z)]  (28) 

which  is  shown  by  proving  that  the  mean  and  NPS  of  each  side  are  equivalent.  The  means  are 
obtained  by  replacing  scatter  with  a  convolution,  and  hence  they  must  be  equal  since  the  convolution 
integral  is  associative.  The  NPS  of  the  left-hand  side  [using  Eq.  (12)]  is 

NPSl(w)  =  [NPS(u)  —  q]  |Ti(u)|2|T2(u)|2  +  q  (29) 

and  that  of  the  right-hand  side  is 

NPSr(u)  =  [NPS(u)  -  q]  |Ti(u)T2(u)|2  +  q.  (30) 

Thus,  they  are  equal  and  Eq.  (28)  is  true  since  |Ti(u)|2|T2(u)|2  =  |Ti(u)T2(u)|2  for  any  complex 
Ti(u)  and  T2(u).  This  means  that  two  cascaded  scattering  processes  can  be  represented  as  a  single 
scattering  process  and  retain  the  same  second-order  statistics  by  using  a  PSF  that  is  a  convolution 
of  the  two  individual  PSFs.  In  addition,  the  order  of  the  two  scattering  processes  in  Eq.  (28)  can  be 
reversed  without  affecting  the  outcome,  showing  that  scatter  commutes  with  itself. 
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6.  The  total  “energy”  in  q(x)  is  given  by  W  as 

/oo 

\q(x)\2dx  =  N  (31) 

-OO 

where  N  is  the  total  number  of  quanta.  This  result  is  equal  also  to  the  total  energy  in  qs{x)  following 
a  scatter  and  hence  scatter  conserves  energy. 


3.  Scatter  in  the  Fourier  Domain 

Scatter  relocates  individual  quanta  by  a  random  vector  A r.  Application  of  the  shift  theorem3  indicates 
that  the  Fourier  transform  of  this  operation  therefore  corresponds  to  multiplication  of  the  Fourier  transform 
of  individual  quanta  by  the  factor  e~i2*k'Ar  which  introduces  a  random  phase  change  and  where  *  represents 
a  vector  dot-product  operation  ( k  *  A r  =  uAx  +  vA y).  If  the  Fourier  transform  of  q(r)  is  Q(k),  where 

Q(k)  =  F  {<z(r)}  =  F  8(r  —  r<)}  =  £  e~i2*kri,  (32) 

L=1  J  »=1 

then  the  Fourier  transform  after  scatter  is  Qs(k)  where 

Qs(k)  =  YJe~i2nkiri+Ari)  (33) 

2—1 

and  where  Arj  is  a  random  vector  variable  describing  the  displacement  of  the  zth  quantum.  It  is  clear  that 
unlike  the  convolution  integral,  the  effect  of  scatter  is  not  multiplicative  in  the  spatial-frequency  domain. 


D.  Relationship  to  the  Detective  Quantum  Efficiency 

The  detective  quantum  efficiency  (DQE)  describes  the  transfer  through  an  imaging  system  of  the  image 
NPS  weighted  by  the  squared  system  transfer  function,  and  can  be  expressed  as28 


DQE(fc)  = 


G2|T(fc)|2NPS(fc) 
NPS  out{k) 


(34) 


where  G  is  the  average  overall  system  gain,  |T(fc)|  is  the  system  MTF,  NPS(fc)  is  the  input  NPS  which  is 
equal  to  q  (the  average  incident  quantum  distribution)  if  the  input  quanta  are  uncorrelated,8  and  NPS0Ut(fe) 
is  the  image  NPS.  Using  this  approach,  the  effect  of  scatter  on  the  DQE  is  obtained  by  setting  G  =  1  and 
noting  that  the  numerator  is  then  equal  to  the  NPS  for  deterministic  blur  NPS d(k)  [Eq.  (11)].  Hence,  the 
effect  of  scatter  on  the  DQE  can  be  expressed  as  the  term  DQE'(fc),  defined  as  the  ratio  of  NPS d{k)  to 
the  actual  NPS,  NPS0Uj(fc): 


DQE'(fc)  = 


NPSd(fc) 
NPS  out{k) 


(35) 


A  deterministic  linear-filter  blur  therefore  has  no  effect  on  the  system  DQE.  The  effect  of  scatter  on  the 
DQE,  DQEj(fc),  is 


DQE's(fc)  = 


NPS  djk) 

NPS  s(k) 

NPSrf(fc)  1 

NPSd(fc)  +  q[  1  -  |T(fc)|2]  g[l  ~  |T(fe)|2] 

NPS(fc)|T(fc)|2 


(36) 

(37) 
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and  is  always  less  than  unity.  For  the  special  case  of  Poisson-distributed  incident  quanta,  NPS(fc)  =  q , 
and  DQE's(fc)  becomes 

DQE.W  =  1  =  |T(fe)|2-  (38) 

1+  |T(fc)P 

Thus,  blur  described  with  a  convolution  does  not  result  in  a  degradation  of  the  DQE,  while  blur  due  to 
scatter  does.  This  result  must  be  interpreted  with  caution  as  in  general  the  DQE  of  a  cascaded  system  is 
not  simply  a  product  of  DQE  factors.29 

Using  this  type  of  linear-systems  approach,  a  more  general  form  of  Eq.  (38)  has  been  derived  elsewhere,29 
giving  the  DQE  of  a  system  consisting  of  an  arbitrary  cascade  of  M  scattering  and  gain  processes  as 


DQE(fc)  = 


1 


M 


1+E 


1  "b  eSi|Ti(fc)|" 


-;Pi(k)=l[gj\Tj(k)\‘2 


j= 1 


2—1 


p  i(k) 


(39) 


where  each  stage  may  be  either  a  scatter  or  a  gain  stage,  but  not  both;  gi  is  the  gain  of  the  ith  stage  having 
a  gain  variance  of  o2.  which  can  also  be  expressed  in  terms  of  a  gain  Poisson  excess,  e9i  —  (<7g  /gi  —  1);  and 
Pj(fc)  is  the  product  of  the  gains  and  squared  MTFs  of  all  stages  up  to  and  including  the  *th  stage.  For 
scattering  stages,  &  =  1,  o2.  =0,  and  e9i  =  -1.  For  gain  stages,  T;(fc)  =  1.  Equation  (39)  describes  the  DQE 
following  a  cascade  of  stochastic  gains  and  scatters  only.  Linear-filter  blurs,  such  as  the  stage  representing 
integration  of  quanta  in  the  detector  elements  as  a  convolution  with  U(x/ax,y/ay),  would  not  normally  be 
included  in  Eq.  (39)  since  the  convolution  integral  does  not  affect  the  DQE. 


CONVOLUTION 


SCATTER 


Notation: 

qd(r)  =  q{r)  *  *  p(r) 

Definition: 

roc 

qd(r)  =  /  q(r')p(r  -  r')  d V 

J  —00 

Fourier  Transform: 

F{qd(r)}  =  F{q(r)}T(k) 

Commutative: 

terms;  self; 

addition;  multiplication 

Distributive: 

addition;  multiplication 

Associative: 

Yes 

Expectation  Value: 

F{qd{r)}  =E{g(r)}**p(r) 

MTF: 

MTF(fc)  = 

NPS: 

NPSd(fe)  =NPS(fe)|T(fc)|2 

Effect  on  DQE: 

none 

Qs{r)  =  q{r )  *s  *5  p(r) 

N 

Qs(r)  =  E^(r  _  tri  + 

2=1 

F{«,(r))=i;e-i2*(r‘+ar‘) 

2=1 

self;  binomial  selection 

No 

No 

E{?s(r)}  =  E{g(r)}**p(r) 

=  E{c/d(r)} 

MTF(fc)  =  |T(fc)| 

NPSs(fc)  =  [NPS(fc)  -  q]\T(k)\2  +  q 
=  NPSd(fc)  +  q\l  —  |T(fc)|2] 
(Rabbani  et  al.) 

system  dependent _ 


Table  1:  Summary  of  the  properties  of  the  (deterministic)  convolution  (*)  and  scatter  (*5)  operators  for 
an  input  image  q(r)  =  <5(r  —  r^),  where  T(fe)  is  the  Fourier  transform  of  p(r). 
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Figure  5:  Simulated  images  showing:  a)  initial  “pre-blurring”  image;  b)  blurred  image  using  the  convolution 
integral;  and  c)  blurred  image  using  the  scatter  operator  with  the  same  PSF  as  b). 

These  characteristics  of  the  scatter  operator  are  summarized  and  compared  with  those  of  the  convolution 
integral  in  Table  I. 


IV.  RESULTS 


A.  Simulated  Images 

The  visual  effect  of  scatter  is  illustrated  by  simulating  images  of  sinusoidal  patterns  of  various  frequencies30 
where  image  blur  is  described  first  as  a  convolution  and  second  as  a  scatter.  It  is  not  intended  that  this 
calculation  describe  a  particular  physical  situation;  rather,  it  is  presented  as  a  visual  comparison. 

A  “pre-blurring”  image  was  generated  assuming  an  average  of  100  detected  quanta  per  0.025-  x  0.025- 
mm2  pixel.  An  isotropic  Gaussian-shaped  PSF,  p(r),  with  a  standard  deviation  width  o  =  4  pixels  (0.1 
mm)  was  used  to  represent  the  blur,  where 


p(r)  = 


1  -r2/2<T2 

2ito2 


(40) 


and  where  r2  =  x2  +  y2.  The  resulting  MTF  associated  with  this  blur  is 


|T(fc)  |  =  e 


-2  TrVfc2 


(41) 


Figure  (5a)  is  the  simulated  pre-blurring  image  on  a  256-  x  512-pixel  matrix.  All  10  sinusoidal  patterns 
can  be  seen,  although  the  image  noise  obscures  some  image  detail.  Figure  (5b)  is  the  result  of  a  conven¬ 
tional  convolution  of  the  pre-blurring  image  with  the  PSF  described  by  Eq.  (40).  It  is  shown  that  only 
approximately  7  patterns  can  be  seen,  and  that  image  noise  has  also  been  reduced  significantly.  Figure  (5c) 
is  the  result  of  a  “stochastic  convolution”  of  the  pre-blurring  image  with  the  same  PSF.  Approximately 
7  patterns  can  still  be  seen,  similar  to  Fig.  (5b),  although  the  noise  in  Fig.  (5c)  is  significantly  greater. 
Theoretically  one  would  expect  both  Figs.  (5b)  and  (5c)  to  have  the  same  spatial  resolution  (same  MTF) 
but  very  different  noise  properties.  The  difference  between  Figs.  (5b)  and  (5c)  illustrates  the  difference 
between  convolution  and  scatter  for  an  input  image  with  Poisson-distributed  quanta. 

The  NPS  measured  from  each  image  in  Fig.  (5)  are  shown  in  Fig.  (6).  Subsequent  to  convolution  with 
the  PSF,  the  NPS  decreases  rapidly  with  increasing  frequency,  and  shows  excellent  agreement  with  the 
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Figure  6:  The  NPS  measured  from  the  image  data  used  in  Fig.  (5)  (thin  lines)  are  compared  with  theoretical 
predictions  (thick  lines)  for  convolution  (deterministic  blur)  and  scatter  operations. 


theoretical  NPS  predicted  by  Eq.  (11).  Subsequent  to  scatter,  the  NPS  is  uniform  in  frequencies,  and 
shows  excellent  agreement  with  the  theoretical  NPS  predicted  by  Eq.  (12). 

B.  Scatter  in  a  Cascaded-Systems  Analysis 

By  extending  the  linear-systems  approach  to  make  use  of  scatter  when  appropriate,  the  linear-systems 
method  can  also  be  used  to  describe  noise  properties  of  these  systems.  In  this  section,  illustrative  examples 
of  this  approach  are  described. 

1.  Image  Focal-Spot  Blur 

The  effect  of  x-ray  focal-spot  blur  on  image  sharpness  has  been  described  by  Rossmann  and  Doi, 31,32 
Burgess,33,34  Barrett  and  Swindell1  and  others  in  terms  of  a  convolution  integral  with  the  appropriate 
focal-spot  PSF.  This  calculation  is  summarized  in  Appendix  II  where  it  is  shown  that 

E{<z(r)}  «  E{4(r)}  *  *  p f(r)  (42) 

where  E{g(r)}  and  E{<?(r)}  represent  the  expectation  (noise-free)  pre-  and  post-blurring  images  respec¬ 
tively,  and  p y(r)  is  the  normalized  focal-spot  PSF  projected  onto  the  image  plane.  If  we  instead  omit 
the  expectation  operator  and  represent  the  pre-blurring  image  as  q{v),  an  actual  sample  distribution  of 
quanta,  then  second-order  statistics  of  the  blurred  image  q{r )  are  preserved  if  we  write 

?(r)«?(r)  *s*sp/(r).  (43) 

The  distribution  q(r)  is  Poisson  distributed  and  hence  the  DQE  associated  with  the  finite  size  of  the 
source  is  given  by  Eq.  (38)  as 

DQE  f(k)  =  |T/(fc)|2  (44) 

where  |T/(fc)|  is  the  focal-spot  MTF.  An  additional  consequence  is  that  while  image  sharpness  is  degraded, 
image  noise  (NPS)  is  unaffected  by  focal-spot  size.  Although  this  result  has  been  known  for  some  time 
(e.g.  Sandrik  and  Wagner13),  it  is  presented  as  an  illustrative  example  to  show  that  this  result  can  also  be 
obtained  using  a  linear-systems  approach  with  the  scatter  operator. 

Off-focal  radiation  is  also  a  source  of  image  quality  degradation.  If  the  source  PSF  is  calculated  to 
include  the  distribution  of  any  off-focal  radiation  contributing  to  image  formation,  the  NPS  and  DQE 
obtained  with  this  approach  will  describe  the  corresponding  image  quality  degradation. 
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2.  Hypothetical  Radiographic  System 

The  hypothetical  system  illustrated  in  Fig.  (1)  is  now  analyzed  using  the  extended  linear-systems 
approach.  Screen  blur  is  represented  with  the  scatter  operator  while  integration  of  light  quanta  in  the 
detector  elements  remains  a  convolution  integral.  The  presampling  detector  signal  therefore  becomes 


d{x,  y)  =  g 


[qi{x,  y)m  *s  *s  pr{x,  y)]P**H 


z,*-) 

0>x  & y  J 


(45) 


This  expression  cannot  be  rearranged  significantly  because  scatter  does  not  commute  with  multiplication 
by  m  (Sec.  III.C.2.ii). 

The  corresponding  NPS  is  obtained  by  cascading  the  NPS  of  each  individual  step  as  summarized  in 
Table  II,  where  the  signal  and  NPS  following  each  stage  are  determined  by  cascading  appropriate  combi¬ 
nations  of  Eqs.  (12)  and  (18).  The  DQE  of  the  system  can  then  be  obtained  using  Eq.  (34).  Equivalently, 
Eq.  (39)  can  be  used  directly  to  obtain 

DQEs(fc)  =  l  +  em  1  -  |Tr(fc)j2  1  -T~ 

+  m  m|Tr(fc)|2  m/3\Tr(k)\2 

(46) 

=  - - - T -  (47) 

i  +  hn  + _ _ _ 

m  mf3\Tr(k)\2 


which  is  the  accepted  result.16,35 

It  should  be  noted  that  if  the  convolution  integral  had  been  (incorrectly)  used  to  describe  image  blur, 
the  DQE  would  have  been 

DQE  d(k)  = - — - —  (48) 


1  +  —  -I — — 

m  mfi 


which  is  a  similar  result  but  lacks  any  spatial-frequency  dependence.  The  difference  between  scatter  and 
deterministic  blur  is  fully  responsible  for  the  frequency  dependence  of  the  DQE. 


STAGE 

IMAGE  QUANTA 

NPS 

X  Rays 

Interacting 

Qr(r) 

NPSj(fc)  =  qi 

Optical  Photons 
Generated 

Qm(r)  =  qi(r)m 

NPSm(fc)  =  m2NPS/(fc)  +  (J^qi 
=  qifh2{l  +  ^) 

Optical  Photons 
Scattered 

Qs(r)  =  qm{r)  *s  *s  Pr(r) 

NPS s(k)  =  [NPSm(fc)  -  gm]|Tr(fc)|2  +  qm 
=  g/m2(  1  +  ^  -  s)|Tr(fc)|2  +  qifh 

Optical  Photons 
Detected 

Qdetir)  =  qs(r)P 

NPS det(k)  =  /32NPSs(fc)  +  P(1  -  P)qs 
=  qim2(  1  +  -  ±)t32\Tr(k)\2  +  qifhp 

Presampling 
Detector  Signal 

ds{x,y)  = 

g0qdet(x,y )  *  *n(^>  %;) 

NPSs(it,  v)  = 

NPSde( (u,  v)gl a2 arsine2 (7rnax)sinc2 (nvay) 

Table  2:  Expressions  describing  the  distribution  of  image  quanta  and  the  NPS  at  each  stage  of  the  hypo¬ 
thetical  system. 
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In  the  absence  of  additive  noise,  the  DQE  of  the  entire  system  including  the  effect  of  focal-spot  blur 
and  screen  quantum  efficiency  can  be  obtained  by  extending  the  model  with  additional  stages: 


d(x,y )  =  gQ 


{< q(x,y )  *s  *5  p  f{r)}am  *s  *s  p  r{x,y)  P  *  *  II  | 


v. 


(49) 


where  q(x,  y)  represents  the  sample  distribution  of  x-ray  quanta  that  would  have  been  obtained  had  there 
been  no  focal-spot  blur  (Appendix  II).  The  DQE  obtained  using  Eq.  (39)  is 


DQE  s(k) 

\  1  —  |T/(fe)|2  l- a  1  +  em  l-|Tr(fc)|2  1-0  1  1 

“  1+  |T/(fe)|2  +  a\Tf(k)\2  +  am\Tf(k)\2+  am\Tf(k)Tr(k)\2  amP\Tf(k)Tr(k)\\ 

=  «IT/(*)12  (50) 

i  ,  fm  j _ \ _ 

m  m/3\Tr(k)\2 

which  includes  the  effect  of  focal-spot  blur. 


V.  DISCUSSION 

The  scatter  operator  described  in  this  article  is  a  translated  point-process,10  and  represents  the  phys¬ 
ical  mechanism  giving  rise  to  the  Rabbani,  Shaw  and  Van  Metter  NPS.  It  is  described  as  a  stochastic 
convolution”  in  which  the  blur  PSF  is  “folded”  into  the  image  signal  through  a  stochastic  mechanism.  It 
is  used  as  an  alternative  to  the  convolution  integral  (linear  filter)  description  of  blur  in  a  linear-systems 
analysis.  Convolution  describes  a  deterministic  cascading  of  weights,  while  the  scatter  operator  describes 
a  cascading  of  probabilities. 

A  deterministic  system  is  one  in  which  the  output  depends  only  on  the  input  (Papoulis  )  and  where 
Eq.  (11)  describes  the  NPS  transfer  through  a  blurring  mechanism.  In  contrast  to  this,  a  stochastic  system 
is  one  in  which  the  system  transfer  function  contains  a  stochastic  element.11  In  the  context  of  this  article, 
this  is  due  to  the  statistical  nature  of  quanta  transferring  image  information  through  the  system.  Physical 
imaging  systems  are  almost  universally  stochastic  systems,  and  thus  the  Rabbani  results  [Eqs.  (12)  and 
(18)]  should  generally  be  used  to  describe  the  NPS  transfer  through  blur  and  gain  processes  in  these 
systems. 

In  the  simplistic  hypothetical  system  discussed  here,  it  was  assumed  that  none  of  the  optical  quanta 
were  absorbed  in  the  screen.  When  this  is  not  a  good  assumption,  an  additional  stage  must  be  introduced 
to  describe  the  fraction  of  generated  optical  quanta  that  are  emitted.  This  is  represented  as  a  binomial 
selection  stage  (a  gain  less  than  unity).  Since  the  scatter  operator  commutes  with  binomial  selection 
(Sec.  III.C.2),  the  relative  order  of  this  selection  and  the  scatter  stage  is  unimportant. 

A  deterministic  blurring  process  has  no  effect  on  the  DQE  of  a  system.  That  is,  both  the  square  of  the 
signal  and  the  NPS  are  passed  through  the  square  of  the  blur  MTF.  This  is  not  the  case  for  a  scatter,  where 
the  DQE  decreases  according  to  Eq.  (12)  and  hence  it  can  be  concluded  that  the  frequency-dependence  of 
the  DQE  is  a  direct  consequence  of  the  scatter  of  image  quanta.  In  the  absence  of  additive  noise,  these 
statistical  properties  -  as  described  by  the  Rabbani  NPS  and  hence  the  scatter  operator  -  fully  account  for 
this  frequency  dependence.  This  observation  would  suggest  that  at  sufficiently  high  spatial  frequencies, 
where  the  MTF  of  a  system  is  necessarily  much  less  than  unity,  scatter  is  always  the  dominant  factor 
responsible  for  degradation  of  the  DQE. 

The  impact  on  the  DQE  of  scatter  in  a  cascaded  system  can  only  be  determined  when  expressed  in  the 
context  of  the  complete  system.  For  instance,  in  the  hypothetical  system  the  DQE  [Eq.  (50)]  is  degraded 
by  scatter  only  when  both  the  terms  |Tr(fc)|2  and  m/?|Tr(fc)|2  are  less  than  unity.  The  second  term 
describes  the  effective  number  of  secondary  optical  quanta  at  a  particular  frequency  that  are  detected  by 
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the  detector  array.  A  value  less  than  unity  indicates  the  presence  of  a  secondary  “quantum  sink”.  Thus, 
it  is  the  secondary  quantum  sink  that  is  responsible  for  making  scatter  a  significant  degrading  factor  in 
the  DQE.  If  the  squared  MTF  term  is  not  small  also,  then  the  secondary  quantum  sink  exists  due  to 
inadequate  gain  rather  than  to  scatter.  A  more  general  comment  can  therefore  be  made  that  for  any 
specified  frequency,  a  scattering  process  degrades  the  system  DQE  significantly  if  the  squared  MTF  of  that 
process  is  significantly  less  than  unity,  and  the  system  has  a  subsequent  secondary  quantum  sink  at  any 
specified  frequency  [Pfik)  <  1  in  Eq.  (39)]. 

VI.  CONCLUSION 

The  convolution  operator  often  used  to  describe  image  blurring  processes  in  linear-systems  theory  passes 
first-order  statistics  (mean  values)  but  not  second-order  statistics  when  applied  to  imaging  systems.  As  a 
consequence,  image  noise  is  generally  not  analyzed  using  linear-systems  theory,  but  rather  with  stochastic 
process  theory.  A  scatter  operator  is  described  here  which  represents  the  physical  process  giving  rise  to 
the  Rabbani  NPS,  and  passes  first-  and  second-order  statistics  thereby  unifying  image  blur  and  image 
noise  calculations  within  a  linear-systems  framework.  A  convolution  blurs  the  input  in  a  deterministic  way 
according  to  the  PSF,  while  scatter  relocates  quanta  according  to  a  probability  given  by  the  same  PSF. 
The  scatter  operator  is  expressed  in  terms  of  random  variables  and  is  consistent  with  a  linear-systems 
approach  to  modeling  linear  and  shift-invariant  imaging  systems,  allowing  for  the  use  of  Fourier-domain 
mathematics. 

Several  properties  of  the  scatter  operator  are  described.  It  is  shown  that  it  commutes  with  itself  and 
with  a  binomial  selection  process  (representing  a  quantum  efficiency).  Strictly  speaking  it  is  not  associative, 
although  a  similar  property  that  is  true  is  described.  The  input  and  output  must  necessarily  be  expressed 
in  units  of  quanta  per  unit  area,  and  the  associated  PSF  must  necessarily  be  normalized  to  unity  area. 
It  is  shown  that  for  LSI  systems  in  general  (in  the  absence  of  additive  noise),  the  transfer  of  second- 
order  statistics  through  scattering  processes  is  fully  responsible  for  the  spatial-frequency  dependence  of  the 
system  DQE. 
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In  this  appendix  it  is  shown  that  the  operation  of  integrating  image  quanta  in  discrete  detector  elements 
of  a  detector  array  can  be  described  with  a  convolution  operation.  The  resulting  presampling  detector 
signal  includes  the  smoothing  effect  caused  by  the  finite  width  of  detector-element  apertures,  but  does  not 
include  aliasing  artifacts.  This  approach  of  using  a  presampling  signal  is  equivalent  to  the  presampling 
optical  transfer  function  described  by  Sones  and  Barnes36  and  Giger  and  Doi24  among  others. 

It  is  assumed  that  the  distribution  of  quanta  that  interact  in  the  detector  can  be  represented  as  q(x ,  y), 
a  generalized  function  consisting  of  the  superposition  of  J-functions  where  each  J-function  represents  a 
single  interacting  quantum.  Thus,  the  signal  from  a  rectangular  detector  element  centered  at  position 
2n,yn  and  having  dimensions  ax  x  ay  is  given  by 
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(A-l  ) 


where  g0  relates  the  detector  signal  (in  detector-signal  units)  to  the  number  of  interacting  quanta. .  Inspec¬ 
tion  of  Eq.  (A-l  )  shows  that  this  integral  can  also  be  expressed  as  the  integral  of  q{x,y)  multiplied  by  a 
rectangular  function  having  unity  height  and  dimensions  ax  and  ay  that  describe  the  element  aperture. 
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Equation  (A-2  )  is  recognized  as  being  the  two-dimensional  cross-correlation  integral  of  q(x,y)  with 
U{x/ax,y/ay),  evaluated  at  position  xn,yn.  Therefore,  we  can  express  the  signal  from  the  nth  detector 

element  as 

dn  =  d{x,y)\x,y=Xn>yn  (A*3  ) 

where  d(x,  y)  is  the  presampling  detector  signal  that  provides  the  detector  element  values  for  all  possible 
element  positions.  When  evaluated  (“sampled”)  at  positions  corresponding  to  the  centers  of  the  actual 
elements,  it  provides  the  set  of  actual  detector-element  values.  The  function  d(x ,  y)  is  given  by 


*«-'«■*•')***{*%)  <a-4) 


in  terms  of  the  two-dimensional  cross-correlation  operator  **,  or  as 

d{x,y)  =9o  q(x,y)**  n  (A_5  ) 

in  terms  of  the  two-dimensional  convolution3  operator  **.  The  aperture  profile  is  symmetric  and  hence 
the  negative  signs  in  Eq.  (A-5  )  can  be  omitted. 


APPENDIX  II:  Focal  Spot  Blur  as  a  Convolution 

In  this  appendix,  a  calculation  showing  that  the  average  effect  of  focal-spot  blur  can  be  expressed  as  a 
convolution  is  summarized  following  Barrett  and  Swindell.1  In  Fig.  (7),  r"  is  a  two-dimensional  coordinate 
in  the  source  plane,  r'  is  in  the  object  plane,  and  r  is  in  the  image  plane  (r"  and  r  have  been  exchanged 
relative  to  Barrett  and  Swindell’s  use.  Care  must  be  taken  because  r,  r'  and  r"  lie  in  three  parallel  but 
different  planes.  If  s(r")  is  an  emission  function  describing  the  emission  of  x  rays  (photons  emitted  per 
unit  area)  from  the  source,  s(r")dV  is  the  mean  number  of  x-ray  photons  emitted  into  all  space  from 
an  elemental  area  d2r".  If  the  x  rays  are  emitted  uniformly  in  all  directions  from  the  source  and  t(r  )  is 
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source,  f(r") 


Figure  7:  Calculation  of  the  focal  spot  MTF. 


the  transmission  factor  through  a  planar  object  at  position  r',  then  the  expectation  x-ray  photon  density 
reaching  the  image  at  position  r  can  be  written  as1 

E{^(r)}  =  -TTa  f  s(r")  cos3 {9)t(r')d2r"  (A-6  ) 

4? Tlyz  J source 

where  6  is  the  angle  of  the  x-ray  beam  with  respect  to  the  source-image  direction,  and  D  is  the  source-image 
distance.  Barrett  and  Swindell  then  show  that  if  D  is  large  with  respect  to  the  source  dimensions, 
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where  m  —  D/d  is  the  geometric  magnification  factor  and  d  is  the  source-object  distance.  It  is  more 
convenient  to  work  with  the  projection  of  s  and  t  onto  the  image  plane.  Therefore,  t([r  —  r0]/m)  is 
replaced  with  ij(r  —  r0)  where  the  subscript  i  indicates  projection  onto  the  image  plane.  The  integral  of 
s(r")  over  all  space  gives  E{IV},  the  expectation  total  number  of  x  rays  emitted  from  the  source.  As  a 
result,  [1  -  m]-2s(r0/[l  -  m])  is  replaced  with  E{7V}P/(r0)  where  P/(r0)  is  the  focal-spot  PSF  projected 
on  to  the  image  plane  (normalized  to  unity  area).  Therefore, 


E{<?(r)}  ~  /_%/M*i(r-r0)d2r0  (A-8  ) 

~  23p/<->  *  * ti^  (a_9  ) 

This  result  shows  that  the  average  blurring  effect  of  the  focal  spot  can  be  represented  as  a  convolution  of 
the  source  profile  with  the  object  transmission  function.  The  commutative  property  of  the  conventional 
convolution  operator  means  that  Eq.  (A-9  )  can  also  be  written  as 

E{q(r)}  fa  ^^(r)  *  *  p/(r)  (A-10  ) 

«  E{(j(r)}  *  *  pf{r).  (A-ll  ) 

where  E{<j(r)}  =  E{Ar}/(47rD2)tl(r)  represents  an  expectation  image  that  has  not  been  blurred  by  the 
focal  spot.  The  expectation  blurred  image  is  therefore  represented  as  the  expectation  “pre-blurrmg”  image 
convolved  with  the  focal-spot  PSF. 
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